c fh Pi Office de la propriete Canadian 
V Li intellectuelle 



du Canada 

Un organisme 
d'lndustrie Canada 



Intellectual Property 
Office 

An Agency of 
Industry Canada 



J * SEPTEMBER 200J H.OQ.O* 



RECd 25 SEP 2003 



WlK) 



PCT 



Certification ./ 
La pr&ente atteste t|ue les documents^'/'. /f'/:;-^-^ 

ci-joints, dont la liste'|^ ; 'at^pHS^hereto and identified below are 

sont des copies audier^i^u^^^dpc^; ^ V ^f^ on file in 

ments d<5pos& au Bureau des brevets. ' the. featent'Qffice. 



Certification 
Tf&s is^o certify that the documents 



Specification and : Dt$wmgi ;as o^^fapa^i^^^^giic^on for Patent Serial No: 
'2,397^89, tKjta^ POULIN, MAREE- 

FLAVBE AUCLAIR^f0RT^RA^ for "Image Model Based on 

N-Pixels and Defined >iti ^gebraib:T6poi0gy/and Applications Thereof. 



1 'N 



PRIORITY DOCUMENT 

SUBMITTED OR TRANSMITTED IN 
COMPLIANCE WITH 
RULE 17.1(a) OR (b) 




CA 02397389 2002-08-09 



-107- 

ABSTRACT OF THE DISCLOSURE 

A data structure in which image is specified by its support, its 
quantities that are linked to the support and allowable generic operations is 
described herein. The algebraic structures of the support are defined using 
algebraic topology concepts such as chain and boundary. The quantities may 
be a scalar, vector, tensor or any other type. They are specified by the cochain. 
The generic operation corresponds to the coboundary operator. The image 
model has several advantages: it allows the derivation of efficient algorithms 
that Operate in any dimension and the unification of mathematics and physics 
to solve classical problems in image processing and computer vision. It can be 
used for binary images as well as for several image acquisition systems. 



CA 02397389 2002-08-09 



-1- 

TITLE OF THE INVENTION 

Image Model based on n-pixels and defined in algebraic 
topology, and applications thereof 

FIELD OF THE INVENTION 

[0001] The present invention relates to an image model based on n- 

pixels. More specifically, the present invention is concerned with a) an image 
model based on n-pixels; b) an application the image model for the resolution 
of diffusion and optical flow; and c) an application of the image model for the 
deformation of curves. 

BACKGROUND OF THE INVENTION 

[0002] People have a notion of what an image is. For instance, for 

psychologists the image is linked to the shape of objects, their depth, the 
relationship of these shapes and their perceptual organization. 

[0003] Artists are focused on how features such as shape, color, 

and perspective are organized to represent a scene that may originate in their 
imagination. 

[0004] Physicists are concerned with the physical phenomena 

produced by a given scene and how they are represented in the image. 

[0005] Neurophysiologists regard images through visual phenomena 

in humans and animals, such as contrast sensitivity, Mach bands, contrast, 
constancy, and depth perception. 

[0006] While a precise definition of the image is elusive, it is clear 

that for certain people an image is the visualization of physiological, perceptual, 
or physical phenomena and for others it is related to semantic content 
Whatever the term image means, the intention is to establish a foundation upon 
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which images of all forms and contents can be discussed with minima! 
confusion. 

[0007] In image processing and computer vision, the foundation is 

related to the understanding of the image formation process. Light generated 
by a source is modified by the objects of the scene. The modified light is 
captured by a system acquisition device, transformed into an appropriate form 
and displayed on a physical support. The content of the resulting image and, 
consequently, its further use, both depend on the properties of the light 
(structured or not, spectral properties such as the range of wavelengths, the 
number of ranges and the intensity), the characteristics of the acquisition 
device, the transformation (analog/digital, pre-processing), the display format 
(temporal or spatial organization of image elements, film, vector, raster). From 
the automatic processing point of view, the image is enhanced to improve its 
perceptual quality or to make some of its features explicit. Usually, its content is 
analyzed via a successive reduction process to construct a more descriptive 
representation in terms of relevant features, which can be used more effectively 
by a decision system (car, robot, etc.) or to help human beings in their daily 
tasks. 

[0008] One of our concern here is to focus on the data structure for 

images and its consequences on the processing scheme. 

[0009] Algebraic topology concepts such as chain, boundary, 

cochain and coboundary are a key to representing images. Algebraic topology 
is well-known domain in mathematics [10, 6, 9]. The literature of algebraic 
topology offer various knowledge that can be used for images. The specialists 
for algebraic topology, however, have made no effort to implement their 
knowledge into computer vision and image processing. Instead of using 
algebraic topology, the specialists on image processing and computer vision 
have limited their selves to develop solutions based on the topology and on 
discrete geometry [7, 8]. 
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|@01©3 Algebraic topology was first introduced into image processing 

and computer vision by Allili and Ziou for topological feature extraction and 
shape representation in binary images [1, 2, 15]. It is used by Allili, Mischaikow 
and Tannenbaum [3] in medical binary images. Auclair et al. [4] used algebraic 
topology for linear and non-linear isotropic diffusion as well as for optical flow in 
gray level images. P. Poulin et al. [12] used algebraic topology for snakes and 
elastic matching in gray level images. 

[0011] In image processing and computer vision, several image 
models have been accepted and are in recurrent use since several decades. 
These image models integrate both the image support and the local quantities 
associated with this support. The image support is formed by pixels. With each 
pixel is associated a scalar quantity called a gray level, or a vector quantity 
called either color at the perceptual level or multispectral at the signal level. 
Existing models differ primarily in terms of how the image support (definition 
and the organization of pixels) and how values are formulated. The well-known 
image models are the function, the random process, and the ordered set The 
image is a function L x *L y -> G m , where L x = {1,..., N x } and L y = {!,..., N y ) , 
N x *N y is the resolution of the image, G = {1 n}, where n is the maximal 

quantity and m the number of image bands. In the case of a binary image n = 
2, image processing has roots in graph theory, language theory, logic, discrete 
geometry, and so on. If n > 1, usually the image is modeled as a real function 
(analog image where G 9 L x ,L y czW). In this case, function theory, functional 

analysis, catastrophe theory, differential equations, and differential geometry 
are the foundation. An image can also be modelled as a collection of random 
variables X{i,j)\(i,j)eL x xL y . In this case, the probability density function, 

moments, sufficient statistics, time series, and the Markov processes are the 
roots. When the image is modelled as an ordered set, discrete mathematics 
and mathematical morphology are the foundation. 
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{0042} Since the introduction of mathematical morphology in sixteen, 

efforts of researchers in these fields are focused on the use of more and more 
complex mathematical, physical or computer concepts as formalism of specific 
problems (edge detection, image segmentation, optical flow, and deformation) 
without questioning the image model. The definition of new image model can 
leads algorithms that are designed on new basis. An image is a physical or 
mathematical quantity where variables (image support) represent geometrical 
or temporal elements such as points, lines, surfaces, and times. For example, 
the image, as a function L x *L y -*G\ can be defined by both the geometrical 
and topological properties of the domains 4x1, and the topological, 
geometrical and analytical properties of G m . Although existing image models 
have deep roots in mathematics or in physics, the variables, the quantity and 
the association between them are not well-defined. For a given computer vision 
or image processing task, no formal mechanism is given for the integration of 
physical, topological, geometrical properties of objects and their behaviors as a 
part of the image model. Consequently, the resulting computational schemes 
are non-modular and sometimes not easy to reproduce. 

OBJECTS OF THE INVENTION 

[0013] An object of the present invention is therefore to provide an 

image model based on n-pixels. 

[0014] Other objects, advantages and features of the present 

invention will become more apparent upon reading of the following non- 
restrictive description of preferred embodiments .thereof, given by way of 
example only with reference to the accompanying drawings. 

[0015] It is to be noted that various documents are referred to 

herein. These documents are hereby incorporated by reference in their 
entirety. 
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BRIEF DESCRIPTION OF TH E DRAWINGS 
[001 6] In the appended drawings: 

[0017] Figure 1 illustrates a 2-cube; the orientation is given by its 

definition [0,1] X [0,1]; 

[0018] Figure 2 illustrates a) a subdivision that is a cubical complex 

and b) a subdivision that is not a cubical complex; 

[0019] Figure 3 illustrates, in solid lines, a primary cubical complex 

and, in dashed lines, a secondary cubical complex; 

[0020] Figure 4 illustrates a) an original image and b) the resulting 

smoothed image; 

[0021] Figure 5 illustrates a body in 3D space at time t. a) Xi is the 

location of the particle Xi, n is a vector normal to the surface, dS is an 
infinitesimal surface patch and dV is an infinitesimal amount of volume, b) f e is 
an external force applied on dS. p is the mass and b is a body force applied on 
dV. c) q is the heat flow passing through dS and r is the heat produced by dV; 

[0022] Figure 6 illustrates three examples of q-pixels in E 2 (h x l 2 ) a) 

0-pixel : li = {2}, l 2 = {1} b) 1-pixel : d = [2,3], l 2 = {1} c) 2-pixel : li = [2,3], l 2 = 
[1.2]; 

[0023] Figure 7 illustrates an example of the coboundary operation; 

[0024] Figure 8 illlustrates a projection onto (a) the tangential part 

and (b)the normal part of the domain; 

[0025] Figure 9 illustrates examples of cochains for one 3-pixel of 

K p1 . a) cochain T b) cochain H ; 

[0026] Figure 10 illustrates examples of cochains a) Q and b) D for 

one 3-pixel of K 8 ; 
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[0027] Figure 11 illustrates examples of cochains a) T and b) G for 

one 3-pixel of K p ; 

[0028] Figure 12 illustrates three different paths between two points; 

[0029] Figure 13 illustrates computational scheme for the unsteady 

problem with no source; 

[0030] Figure 14 illustrates two cubical complexes for a 5 X 5 image; 

[0031] Figure 1 5 illustrates y p , one 2-pixel of K p ; 

[0032] Figure 16 illustrates a) y E (in dashed lines) one 2-pixel of K 3 

intersecting four pixels of K p ; b) Cochains T and G. c) Cochain Q; 

[0033] Figure 17 illustrates a) One 3-pixel of K tf surrounding one 1- 

pixelofK p '; 

[0034] Figure 1 8 illustrates Xs for one 2-pixel of K p ; 

[0035] Figure 19 illustrates one 3-pixel of K* intersecting four 3- 

pixels of a) Cochain T; b) Cochain G; c) Cochain Q; 

[0036] Figure 20 illustrates two 2-pixels of K s with X's; 

[0037] Figure 21 illustrates a) Physics-based linear isotropic 

diffusion, o = {2.0, 4.0. 5.0}; b) Convolved, for same scales; 

[0038] Figure 22 illustrates first images of three sequences, (a) 

Rotating sphere sequence; (b) Hamburg taxi sequence; and (c) tree sequence; 

[0039] Figure 23 illustrates a flow pattern computed for the rotating 

for the rotating sphere sequence, using a) physics-based method and b) 
iterative finite difference method; 
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[0040] Figure 24 illustrates a flow pattern computed for the Hamburg 

taxi sequence, using a) physics-based method and b) iterative finite difference 
method; 

[0041] Figure 25 illustrates a flow pattern computed for the tree 

sequence, using a) physics-based method and b) iterative finite difference 
method; 

[0042] Figure 26 illustrates a first image of tree sequence with white 

noise added (standard deviation of 10); 

[0043] Figure 27: Flow pattern computed for the tree sequence with 

white noise added, using a) physics-based method and b) iterative finite 
difference method; 

[0044] Figure 28 illustrates a section of the peppers image (o = 5); a) 

Original image with noise added; b) PB method; and c) FD method; 

[0045] Figure 29 illustrates a section of the peppers image (o = 5); a) 

Original image; b) Original with white noise added; c) PB method; d) FD 
method; 

[0046] Figure 30 illustrates: a) PB, o = {1.0, 3.0, 5.0}; b) FD; for 
same scales; 

[0047] Figure 31 illustrates: a) Original image; b) With noise added 
(std dev =10); 

[0048] Figure 32 illustrates: a) PB, o = {4.0, 8.0}; b) FD; for same 
scales; 

[0049] Figure 33 illustrates: a) PB and (b) FD methods with a « 8.0; 



[0050] Figure 34 illustrates a body of arbitrary size, shape and 

material. AS is a surface patch, f is a vector of surface forces applied on AS, 
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AV is an amount of volume with a mass p, b is a vector of body forces applied 
on AV; 

[0051] Figure 35 illustrates a cutting plane passing through a point 

and partitioning the body into two sections; 

[0052] Figure 36 illustrates a force Af acting on a cutting plane with 

normal vector n; 

[0053] Figure 37 shows examples of stresses on a body: (a) 

Original body; (b) Normal stress after an extension of the body; (c) Normal 
stress after a compression of the body; and (d) Shear stress after a distortion of 
the body; 

[0054] Figure 38 illustrates the component of the stress in the 

direction of xr, 

[0055] Figure 39 shows the deformation (a) and distortion (b) of a 

body subjected to stresses; in both figures, the rectangle ABCD has been 
deformed or sheared into A'B'CD; 

[0056] Figure 40 illustrates normal strain of a body; 

[0057] Figure 41 illustrates shear strain in a body; 

[0058] Figure 42: A body B in motion subjected to forces; ti n and pbj 

are respectively the traction and body forces in the direction of xi; 

[0059] Figure 43 illustrates (a) the kinematic equation; (b) the 

constitutive equation; and (c) the conservation equation; 

[0060] Figure 44 illustrates the decomposition of the linear elasticity 

problem into basic laws; 
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[0061] Figure 45 shows three examples of q-pixels in R 2 (h x l 2 ) a) 0- 

pixel : It = {2}, l 2 = {1} b) 1-pixel : h = [2,3], l 2 = {1} c) 2-pixel : li = [2,3], I 2 = 
[1.2]; 

[0062] Figure 46 is an example of the coboundary operation; 

[0063] Figure 47 illustrates examples of cochains for a 3-pixel of K p ; 

(a) Cochain V (b) Cochain D; 

[0064] Figure 48 illustrates examples of cochains S (a) and F (b) for 

a 3-pixel of K p ; 

[0065] Figure 49 shows two cubical complexes for a 5 x 5 image; 

[0066] Figure 50 illustrates a 2-pixel of K p and the topological 

quantities associated with it; 

[0067] Figure 51 illustrates: (a) yf in dashed lines (b) 2-cochains D 

and 14 (c) 2-cochain S; 

[0068] Figure 52 illustrates five adjacent vertices in a curve and its 

deformation; 

[0069] Figure 53 is a table that shows typical values of the Poisson's 

ratio and Young's modulus of elasticity for some materials; 

[0070] Figure 54 shows: (a) Initial curve (b) Bright line plausibility 

image; 

[0071] Figure 55 is the results for an aerial image; (a) PB method 

(b) FEM method; 

[0072] Figure 56 is an initial curve for a SAR image; 

[0073] Figure 57 shows line plausibility for a SAR image; 
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[00741 

PB method; 

[00751 

FEM method; 

[0076] 

[0077J 

[0078] 

method); 

[0079] 
[0080] 



Figure 58 illustrates a road correction for a SAR image with 

Figure 59 illustrates a road correction for a SAR image with 

Figure 60 shows an initial curve; 

Figure 61 shows a bright line plausibility image; 

Figure 62 illustrates corrections for a multiband image (PB 

Figure 63 illustrates corrections for a Landsat 7 image (FEM); 
Figure 64 illustrates initial (a) and corrected (b) curves for a 



synthetic image; and 

[0081] Figure 65 (a) through (c): shape recovery of a curve when the 

external forces are removed, (d) Both initial curve (in white) and final curve (in 
black). 

DESCRIPTION OF THE EMBODIMENT 

[0082] An aspect of the present invention is generally based on a 

data structure for images based on n-pixels, in which the image support, the 
image quantity and the allowable operations are specified separately. In this 
data structure, mathematics and physics are unified; that is, the data structure 
allows taking into account constraints originating in physics, mathematics, and 
the future use of the image. The image dimension is explicit, which allows us to 
design algorithms that operate in any dimension. Other advantages and new 
open problems of this data structure are discussed hereinbelow. 



[0083] One of the goals of the present invention is to give a 

computational image model or a data structure that should be capable of 
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capturing all objects properties that are needed for the given task. A data 
structure of the image is the formal specification of the image variables, image 
quantity, the association between quantities and variables that allows to 
capture the geometrical and topological properties of objects as well as their 
physical and mathematical behavior. The abstract view of image as a data 
structure as is used in computer programs is defined by the attributes and a 
collection of meaningful operations. The attributes are the image support and 
quantities that are assigned to the image support such as the image radiometry 
(e.g., color and grey level) or any feature that can be deduced from the 
radiometry (e.g. texture). These quantities are scalar, vector or tensor. The 
allowable operations are of two kinds: the operations that are problem 
independent such as read and write and those that are problem defendant such 
as object deformation, diffusion and optical flow. To summarize, the image is 
defined by the support, quantities and allowable operations. We will explain in 
details each of these three items hereinbelow. 

Dinmagj© support 

[0084] Often, the discretization of an image is achieved via a cubic 

tessellation. For example, a 2D image support is formed by unit squares 
commonly called pixels. Similarly, a 3D image is a tessellation of unit cubes 
commonly called voxels. More generally, the present invention will consider the 
image in n dimensions as a set of unit n-cubes, which we will name n-pixels. 
When n = 0, the image is set of points; when n = 1, a set of edges; when n = 2, 
a set of squares, and so on. Any two n-pixe!s are either disjoint or intersect 
through a common i-pixel, where i < n. This subdivision of the image support is 
not unique. Several other geometrical forms such as, for example, triangles or 
hexagons can be used. It has been proven that the topological features of the 
image support do not depend on the subdivision used [13]. The cubical 
subdivision is commonly used in image processing and computer vision and 
will therefore be used herein. One does not need to explicitly define an 
orientation on pixels. Indeed, the definition of a pixel as a product of intervals in 
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a certain order coupled with the natural order of real numbers imposes an 
orientation on each coordinate axis and also on the canonical basis of 9T (see 
Figure 1). 

[0085] Formally, a pixel a e« n as a geometric entity is translated 

into an algebraic structure as follows: 

<T = ii X J 2 X X I n 

where x is the Cartesian product and l f is either a singleton or an interval of unit 
length with integer endpoints; i.e., h is either the singleton {*} , in which case it 
is said to be a degenerate interval, or the closed interval [M+l] for some 
integer k. The number q e {0,1,..., n) of non-degenerate intervals in equation 1 is 
by definition the dimension of a , which will be referred to as a q-pixel. For 
qZl, let 7 = {k 0 ,k l ,...,k q . I }be the ordered subset {l,2,...,n}of indices for which 

/ ^[a^bj] is non-degenerate. Define 

A kj a = h x . . . x I kj -i x {a,j} x I k .+i x . . . x I n 

and 

B kj a = /iX.,.x I k .-i x {bj} x I kj+ i x . . . x I n 

A t a and B k a are, respectively, the front (q-1)-face and the back (q-i)-face of 

a . Each of these (q-1}4aces is a (q-1)-pixel. These faces are then called (q-1> 
faces of a . In the same manner, one can define the (q-2)-faces, ... t down to 
the 0-faces of a. In Figure 1, we see a 2-pixel A, with its 1 -faces a, b, c, d. 
The 0-faces of A are the vertices that are not represented here for the sake of 
clarity of the picture. The boundary of the q-pixel a allows us to write the 
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relationship between a q-pixel and its (q-1)-faces in algebraic form. It is by 
definition the alternating sum of its (q-l)-faces, i.e. 

i=o (2) 

[0086] For example, the boundary of A in Figure 1 is 

d 2 a m (-0 x [0,1] + 1 x [0,1]) - ([0,1] x 0 - [0,1] x 1) = (c - a) + (d - b) . 

[0087] So far, the pixel, its faces, and the association between them 

have been defined geometrically and algebraically. We will now define the 
image support as a geometrical entity, called a cubical complex, and then write 
its algebraic structure. A cubical complex K in 9T is a collection of q-pixels 
where 0 <> q < n such that : 

1 . Every face of a pixel in K is also in K. 

2. The intersection of any two pixels of K is either empty or a face of each of 
them. 

[0088] The first condition implies that all faces of a pixel belong to 

the cubical complex. 

[0089] The second condition concerns the organization of the 

cubical complex. If the intersection of two pixels is empty, then the image 
support is formed by one or more connected components. This condition 
provides an image support that is more general than existing definitions since it 
allows a formal specification of a cubical complex that is formed either by one 
or more image supports (e.g., an image sequence) or by several distinct binary 
objects. When the intersection is a face, certain geometrical configurations of 
the complex are ruled out. For example, Figure 2a illustrates a subdivision that 
is a cubical complex and Figure 2b illustrates a subdivision that is not a cubical 
complex. 
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[0090] The dimension of K is by definition the largest number q for 

which K contains a q-pixel. 

[0091] As in the case of the q-pixel, the cubical complex can be 

written in algebraic form. Given a topological space X eST in terms of a 
cubical complex, we denote E q ={<rj,...,<r^}. the set of all q-pixels of X. A q- 
chain in X is a formal sum of integer multiples of elements of E q . More 
precisely, it is a linear combination 

AioJ + .-. + A^of' (3) 

where A, A^are integers. For example, in Figure 1, a-c+d-b is a 1 -chain. 

Two q-chains are added by adding corresponding coefficients. The set of q- 
chains can be given the structure of a free Abelian group with basis E q , usually 
denoted by C 9 (X). Since we will only be concerned with finite complexes, the 
groups C q (X) are finitely generated and C q (X) = 0 if q is greater than the 
dimension of the complex; naturally, C 9 (X) = 0 if q < 0. 

[0092] To define the chains that ate associated with the faces of 

pixels of a cubical complex, we extend equation 2 by linearity to all q-chains 
and obtain the boundary map 5, : C^X)-*^ (X). Note that S 0 =0 since 
C. y (X)~0. The boundary map satisfies the following very important property 
[9]: 

d q O d q+ l = 0 (4) 

[0093] To summarize, the discrete image support of any dimension 

n is formed by n-plxels. Unlike conventional image models, the n-pixel is a 
dimensional geometric entity formed by other geometric entities called faces. 
The geometrical n-pixel is translated into a recurrent algebraic structure, more 
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reliable for mathematical handling. All n-pixels of the Image support form a 
cubical complex, a geometrical entity that is translated into an algebraic 
structure called the chain. The association between the n-pixel and its faces, 
and hence between chains of successive dimensions is given by the boundary 
operator. It is to be noted that the use of any other geometrical primitive such 
as the triangle or hexagon for subdividing the image support affects the 
computational complexity of the derived algorithms, but it affects neither the 
topological features of the image support nor the computation rules for these 
features. 

Image Quantities 

[0094] We introduced hereinabove the concept of the finite cubical 

complex to give an algebraic description of the discrete image support. A 
similar formulation is required to describe the image field (scalar, vector, matrix) 
over the discrete image support. For this purpose, we return to the chain 
concept described above and give a more general definition. Let us consider a 
cubical complex K of dimension n. We associate with each q-pixel (g £ n ) of K 
a coefficient in the ring (A,+.*), where the elements may be scalars (gray level), 
vectors (color, gradient), matrices (Hessian), etc. The chain is the formal sum 
Ytfac', • wnere *t is a coe fficient in (A.+ *). and the generators c\, Vi form a 
basis of an Abelian group. The chain can be seen as a vector 
where iV,is the number of generators used. Consequently, two chains can be 
summed by adding their coefficients and multiplied using the scalar product. 
The addition and multiplication are taken according to the definition of a ring. 
Moreover, one can define the null chain (resp. unit chain) whose coefficients 
are all equal to the null (resp. unit) element of the + (resp. *) operation of 
(A.+,*). Consequently, q-chains define an image model that has attractive 
computational properties since they form a rich algebraic structure, the module 
(i.e. a vector space defined on a ring). 
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[0095] To briefly show how to use the chain model in image 

processing, we present a simple example concerning any global transform of 
the image such as histogram equalization or thresholding. The chain 
coefficients are the gray levels. The formal expression of the global transform 
implies two applications: ff:(f4,+,*) and H : 04,+,*) -> (A 9 + 9 *) , 
where (^,+ *) is a module. They are defined by HC^fcc.) = 2/f(^)q . 

[0096] The drawback of the conventional chain-based image model 

is that the physical or mathematical quantities and the image support are 
described together in the formal sum. Consequently, the chain coefficients 
combine mathematical or physical quantities, pixel orientation, and possibly 
other quantities such as weights associated with pixels. For example, there is 
ambiguity in the interpretation of the sign of A f : it may be due to the orientation, 
the physical quantities, the weights or their multiplication. This image model 
can be considered adequate, especially as in physics [14], engineering and 
computer graphics [11, 5}, chains have been used to model quantities. 
However, we propose to refine this quantity model to overcome the confusion 
produced by the chain coefficients. 

[0097] To reflect only the geometry of the image support (e.g., 

orientation and multiplicity), in what follows we consider n-chains with integer 
coefficients. We are looking for an application F q : C q (X) -> (5,+,*) , which 

associates a global quantity (energy, gray level, color, flux, tensor, etc.) with all 
q-pixels, where q < n and (£>+,*) is a ring. As in the case of the chain-based 

image model, for two adjacent q-pixels c\ and c\ % J 5 ; must satisfy 

^,(Vj +^<%)=^ F q( c \)+^ F q( c q)' which means that the sum of the quantities 
generated within each q-pixel is equal to the quantities generated within the two 
q-pixels. F q can be extended by linearity to any q-chain ]£A C 9 » wh ®re each \ 
is an integer, as follows: 
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[0098] F 9 is called a q-cochain. To illustrate the cochain concept, let 

us consider a 2-pixel ca and a vector field V. A 0-cochain is defined by the 
value of V at 0-pixels. A 1 -cochain is jjv.ds , the line integral over the faces of 

c 2 . A 2-cochain is ^VJS , the surface integral over the 2-plxel c 2 , and the V the 

dot product. 

[0099] To summarize, a cochain allows us to associate quantities 

with an n-pixel and with its faces. Unlike existing image models, our model 
provides a rich structure that allows the definition of both local and global 
quantities. 

Operations 

[00100] A generic operation that can be instantiated depending on the 

problem that we are dealing with will now be defined. Having in mind that the q- 
pixel has 3 Q faces, the generic operation should specify the algebraic 
relationship between the quantities (i.e.. cochains) associated with these faces. 
Based on the linearity principle, the quantity of a given q-pixel is transferred to 
its cofaces with the same or opposite sign, according to the agreement of its 
orientation and the orientation of its cofaces. The quantities that are transferred 
to the (q+1)-pixel by its faces are summed. More formally, it should be recalled 
that the relationship between the (q-1)-chain and the q-chain is given by the 
boundary operator. Similarly, the relationship between the q-cochain and the 
(q+1)-cochain is given by the coboundary operator S q :C 9 ->C* + \ where C q is 
the Abelian group of q-cochains. Given a (q+1)-chain c, this operator is defined 
by: 



8 q F q (c) = F Q (d q+1 c) 



(5) 
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[001 01 J The capacity of the cochain and the coboundary to model a 

given problem will now be discussed. The cochain is a linear application and 
should fulfill the coboundary requirement in equation 5. Thus a question 
concerns the limits in modeling a given quantity. It is difficult to answer this 
question for the general case. Much investigation is needed first We only 
answer by identifying several problems that can be modeled by the cochain 
and coboundary. Certain problems such as convolution and its applications 
(smoothing, numerical differentiation, high-pass filtering, noise estimation) and 
the Fourier transform can be modeled by the cochain without approximation 
since they only require setting the coefficients A^of the cochain to specific 

values (see [16] for the case of numerical differentiation). Thresholding can be 
represented by the cochain without approximation since 
H CILfaQi) = » where H : (5,+,*) -» (5,+,*) . Other problems can be 

broken down to basic laws, each of which can be described by the topological 
equation in (equation 5). For example, it has been pointed out [14] that many 
physics problems can be broken down into basic physical laws such as balance 
and constitutive laws. As we will show hereinbelow, balance law can be written 
in a discrete form by using the topological equation in (equation 5) without 
approximation. The constitutive laws cannot be translated in algebraic form 
without approximation. Usually, they requires the link between cochains that 
belong to different cubical complexes. For example, in the case of dual 
complexes, two cochains are linked by an algebraic linear system. We call this 
transformation the codual operation. More generally, 0-cochains represent 
local quantities. However, q-cochains (q >1) represent global quantities (e.g., 
an integral or the summation of a differential form) since they are associated 
with the algebraic structure of an edge, an area, a volume, etc. Thus, the 
cochain, coboundary, and codual are generic algebraic structures that can be 
instantiated by physical or mathematical laws. The exact translation of given 
problem in terms of q-cochains is possible if one is able to find the basic laws 
that can be described without approximation by either cochains or the 
topological equation in (equation 5). 
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[00102] In our previous example, we can interpret the coboundary as 

follows: assuming that the vector field is conservative F = Av, we have 
|Av.<fe = v(6)-v(fl) which is a coboundary operator since it may be written as 

^o fT o( c i) = J o( s i c i)» where a and b are the faces of c n . Similarly, 
jdiv(V)dS = jv.nds , where n is the outward unit normal vector to dc 2 , which is 

e, dc 2 

a coboundary operator since it may be written as S^(c 2 )^F^d 2 c 2 ). This 
example shows that the coboundary operator may be an exact discrete 
representation of the fundamental calculus theorems (line integral and Gauss). 

[00103] Thanks to the chain, cochain, coboundary and codual 

concepts, the image model of the present invention can take into account the 
mathematical or physical laws related to the application. It can thus be 
instantiated to the various problems of image processing and computer vision. 
The underlying computational framework is strongly different form existing 
frameworks. For example, let us consider physics based problems such as 
optical flow, diffusion and deformation. Existing frameworks can be 
summarized as follow: 1) modeling by partial differential equation; 2) resolving 
the PDE by using numerical analysis scheme or Fourier space. 

[00104] The computational framework, according to an aspect of the 

present invention, can be summarized as follow: 1) identification of its basic 
laws; 2) the definition of the image support that is the number of cubical 
complexes and the dimension of the cubes (e.g., for the case of a multi- 
resolution processing); 3) the definition of global and local quantities; 5) 
instantiation of the coboundary and codual operators; 6) the resolution of 
resulting algebraic system. The advantages of this framework are described 
hereinbelow and will be better defined in two practical examples. 

[00105] To summarize, the coboundary operator links quantities 

associated with the faces of an n-pixel. Codual operator links quantities 
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associated to complexes of an image. If a given problem can be broken down 
into basic laws, the cochain and coboundary are the discrete representation of 
these basic laws. Cochain, coboundary, and codual are generic concepts that 
can be instantiated by physical or mathematical laws. Thus, they can be used 
in various computer vision and image processing problems. 

An Example 

[00106] Let us consider the linear isotropic diffusion in gray level 

images which is a physics-based problem. The reader can find all details in [4]. 
For the sake of clarity and without loss of generality, we will limit ourselves to 
considering the 2D global differential equation for heat flow in a homogeneous 
medium. Let V be a 3D region with boundary S inside the flow. The rate of heat 
flow across S out of V is given by: 

///^(x,t)«v = -/// v V.(Avr ( x, i )Mv (6) 

where x is a spatial vector, t time, and X a positive constant 

[00107] To resolve this problem, the image is defined by a continuous 

scalar field T, the temperature (i.e., gray level), two dual cubical complexes 
(i.e., two chains), and three cochains. If only one cubical complex is used, two 
different orientations are associated with each 1-face. To overcome this 
problem, two dual complexes (primary and secondaiy) are used (see Figure 3). 
Concerning the use of three cochains, it has been pointed out in [4] that this 
EDP is formed by three basic physical laws. Each cochain is associated with 
one law. The first is the thermal tension law (also known as Fick's Law), which 
states that heat flows from regions of higher temperature towards regions of 
lower temperature. The direction of the gradient VT is the direction of the 
largest increase in temperature; the heat flows in the opposite direction. 
Formally, this law is written as follows: 



m 
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ff (x,*) = -VT(x,*) 



(7) 



[00108] The primary cubical complex is the support for this 

balance law. The orientation plotted on this cubical complex is the direction of 
the path on which the integral is computed. Let us assume that the 0-cochain is 
the temperature T associated with 0-pixel Co. The 1-cochain G associated with 
1 -pixel Ci is the global thermal tension transferred by the two 0-pixels that are 
the faces of ci. Consequently, the topological equation is: 



GM = <5oT(d) = T{d x c x ) 



(8) 



[00109] By the linearity of cochains, this topological equation is valid 

for all 1 -chains. 

[00110] The second law, called the heat source law, concerns the net 

outflow of internal thermal energy at the point x and time L This is a balance 
law. It is given by: 



<r(x,*) = V.g(x,i) 



(9) 



[00111] When V.q(x,t) > 0, the outflow is positive and thermal energy 

must flow away from x. Similarly, if V.q(x,t) < 0, the inflow is larger than the 
outflow and thermal energy increases at x. The equilibrium for a diffusion 
process is attained if V.q(x,t) = 0. Let us consider the secondary cubical 
complex. The orientation plotted on this cubical complex is the direction of the 
flow. The 2-cochain Z associated with the 2-pixel c 2 is the global heat 
transferred by the faces of 



£(<*) = 6 l Q(c 2 ) = Q(d 2 c 2 ) 
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[00112] This topological equation is valid for alt 2-chains. The third 

law is constitutive (it depends on the medium feature). It makes the link 
between the flow density law and the heat source law. It is given by: 

q( x ,t) = \g(x,t) 

[00113] This equation cannot be discretized exactly, since it involves 

global quantities defined on two dual complexes. Consequently, the 2-cochain 
± cannot be computed without approximation from the 1 -cochain G. For 
example, they can be linked as a linear equation system A0 = ±, where A is 
the coefficient matrix. Figure 4 gives an original image and the image smoothed 
using this computational scheme. 

[00114] The data structure associated to the linear diffusion problem 

is defined by: 1) two dual cubical complexes; 2) two cochains G and A for 
global quantities and a scalar field T; 3) two coboundary operations for balance 
laws and an codual operation that represent the constitutive law. The 
framework is summarized as follows: g is approximated by a bilinear 
polynomial. The cochain G is computed using a line integral, q is computed 
using the constitutive law in equation 1 1 . The cochain ± is computed using 
Gauss's theorem. Finally, a system of linear equations obtained from equation 
1 1 is obtained where the unknowns are T. It should be noted that the cochains 
in equations 8 and 10 are computed without approximation. 

[00115] The image model according to an aspect of the present 

invention and described hereinabove may generally be characterized by three 
major points: 1) the image support and quantities are defined separately and 
then linked together via algebraic language; 2) the pixel is dimensional and is 
written in an algebraic form; 3) both local and global quantities are represented 
by the cochains and coboundary operator. 
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[00116] Each of these specificities will now be discussed to show 

their straightforward consequences for image processing. 

[00117] The separability of the image model allows the distinction 

between image variables and image quantities. The image variables offers 
numerous possibilities for existing mathematical formulations such as the use 
of algebraic topology to help in the design of algorithms. For example, binary 
image algorithms are written as algebraic systems [16]. The well defined 
quantities allow the use of physics, vector analysis or differential forms in the 
design of algorithms. Taking image support and image quantities together, well- 
known problems such as those of diffusion and optical flow in gray level images 
can be written as algebraic systems. Furthermore, the transfer of quantities 
between a given domain and its boundary is straightforward, using the 
concepts of cochain and coboundary as a general framework. For example, as 
we have shown, in vector calculus, this transfer is easier thanks to the three 
fundamental calculus theorems, the line integral, Stokes's theorem, and 
Gausses theorem, 

[00118] Unlike existing image models, by considering the pixel as 

dimensional primitive the connectivity paradox of the image support is avoided 
[8]. That is, the well-known Jordan theorem is fulfilled. Its decomposition into 
faces and the use of cubical structures such as chains make the dimension of 
the image explicit. Algorithms designed according to this formalism operate in 
any dimension. This fact overcomes the traditional limitations that we face in 
designing an algorithm, say in one dimension, extending it to two dimensions, 
and then to three dimensions, and so on. Each extension step may be a difficult 
task. 

[00119] The definition of the cochain depends on the problem that we 

are dealing with. Thus, this image model offers real flexibility for the integration 
of mathematical objects (scalar, vector, tensor) and physical laws (balance, 
constitutive). Furthermore, the use of global quantities associated with an n- 
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pixel implies noise reduction. In fact, global quantities are computed from the 
field by using the integral or the discrete summation. As the opposite operation 
from the differentiation, which enhances high frequencies of the image, the 
integral performs a smoothing operation. It allows us to reduce the order of the 
derivative used in an image-processing scheme. Consequently, the introduction 
of global quantities may allow the use of higher-order derivatives. 

[00120] Another contribution of the image model according to an 

aspect of the present invention concerns the numerical scheme used to solve 
nonlinear problems such as the diffusion problem and elastic matching. It 
should be recalled that, usually, a problem is formulated and then a numerical 
analysis scheme is used. The numerical analysis scheme may not have been 
derived for exactly this formulation and many approximations must be made. 
The explanation of intermediate results is not available. Consequently, no clear 
idea is available about the convergence of the numerical analysis scheme and 
the numerical results obtained may be a broad approximation of the desired 
solution. Based on the problems tackled, we conclude that in the image model 
presented here, the numerical scheme is deduced from the problem model with 
little or no approximation [4, 12]. In fact, various problems may be broken down 
into basic laws and then reformulated in terms of cochains and coboundaries. 
Thus they can be written as linear algebraic systems and solved. 

PRACTICAL EXAMPLE #1: 

PHYSICS-BASED RESOLUTION OF DIFFUSION AND OPTICAL FLOW 

[00121] We will now describe an alternative to partial differential 

equations (PDEs) for the solution of three problems in computer vision: linear 
isotropic diffusion, optical flow and nonlinear diffusion. These three problems 
are modeled using the heat transfer problem. Traditionally, the method for 
solving physics-based problems such as heat transfer is to discretfze and solve 
a PDE by a purely mathematical process. Instead of using the PDE, we 
propose to decompose the global heat problem into basic laws. We show that 
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some of the basic laws admit ah exact global version since they arise from 
conservation principles. We also show that the assumptions made on the other 
basic laws can be made wisely, taking into account knowledge about the 
problem and the domain. We use the above-described image model which 
allows us to encode physical laws by linking a global value on a domain with 
values on its boundary. The resulting algorithm performs in any dimension. The 
numerical scheme is derived in a straightforward way from the problem 
modeled. It thus provides a physical explanation of each solving step in the 
solution. 

Background 

[00122] In recent years, partial differential equations (PDEs) have 

attracted increasing interest in the field of computer vision. Since PDEs have 
been the subject of much study by numericians, powerful numerical schemes 
have been developed to solve them. Consequently, domains such as image 
enhancement, restoration, multi-scale analysis and surface evolution all benefit 
greatly from PDEs [25]. 

[00123] One important class of equations governing certain physical 

processes is the linear elliptic PDE of the general form known as the Helmholz 
equation: 

VMx) + p(x)u(x) = /(x) (12) 

where x denotes a vector in the n-dimensional space, u(x) is the dependent 
variable, V 2 is the Laplacian operator, and p(x) and f(x) are spatial functions. 
When p(x) = 0, we have the Poisson equation [21 , 32] (also known as the non- 
homogeneous Laplace equation [30, 44]). One of the physical processes 
governed by Equation 12 is steady-state heat transfer. 

[00124] In the field of computer vision, Equation 12 may arise from 

two approaches. The first is variational calculus. As a matter of fact, many 
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problems such as shape from shading [39], surface reconstruction [32 f 40] and 
the computation of optical flow 129] can be formulated as variational problems. 
The solutions to these variational problems are given by Euler-Lagrenge 
equations, which are in the form of Equation 12 [31 , 32]. The second approach 
is physics-based. For example, diffusion processes arise from heat equations 
and shock filters from work in fluid mechanics [25]. For both the variational and 
the physics-based approaches, the resulting PDEs are continuous and have to 
be discretized. 

[00125] Traditionally, the discretization of PDEs in computer vision 

has been done by applying finite difference methods [23, 31, 39, 40], Equation 
12 is solved iteratively using either a direct Fourier-based Poisson solver for 
each iteration [39], finite elements [24], or spectral methods [32], Iterative 
methods such as those in [39] do not ensure convergence unless smoothness 
is very high [21]. 

[00126] We can summarize the existing methods for the resolution of 

problems involving PDEs as follows: 1) identification of basic laws; 2) 
combination of the basic laws in order to write the PDEs; 3) discretization of the 
PDEs; 4) resolution of the PDEs via a numerical method. 

[00127] This process, which has been used in various fields of 

application, is purely mathematical. Consequently, it has the following 
drawbacks: 1) Some quantities involved in the solution process do not have a 
physical interpretation; 2) This lack of interpretation is manifested in 
intermediate solutions involving iterative processes and since these solutions 
cannot be physically explained, discovery of the optimal solution cannot be 
ensured in an optimal time. 



m 
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Solution 

[00128] To overcome these drawbacks, an alternative to PDE 

resolution in the context of the heat transfer problem is proposed and will be 
described hereinbelow. 

[00129] Generally, basic laws in physics-based problems are 

combined into a global conservation equation [42] that is valid on the whole 
body or a part of it. A local conservation equation (PDE) is then obtained by 
considering the particular case of a part of a body reduced to a point 

[00130] In discrete problems such as those encountered in computer 

vision, the continuous domain is subdivided into many sub-domains in which 
there is only one value available, which can be considered as a global value. 
Therefore, instead of using the PDE, we propose to use the global conservation 
equation directly on each sub-domain. 

[00131] In order to handle these physical laws which link global 

values at points, lines, surfaces, volumes, etc. we use the image model with 
roots in computational algebraic topology described hereinabove. This model 
makes it possible to represent global values on entities of any dimension at the 
same time. 

[00132] The above described methodology presents a number of 

advantages: 

1) Many of the basic laws arise out of conservation principles and hence 
they are valid either at a point (local form as in Equation 12) or over an 
entire region (global form). Fundamental theorems of calculus such as 
the Gauss, Green and Line Integral theorems allow the computation of 
the coboundary operator without any approximation. 

2) Some laws require approximations that can be performed wisely, taking 
into account knowledge about the problem and the domain. 
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3) The intermediate results have a physical explanation because they 
represent physical quantities. For that reason, every step has a physical 
interpretation. Thus we no longer have problems of non-optimality of the 
solution, because we avoid non-temporal iterative methods. 

4) As mentioned earlier, this method can be used with other physics based 
problems by applying the appropriate basic laws [36]. 

5) Thanks to the image model, the resulting algorithm performs in any 
dimension. 

6) The computational rules associated with the coboundary operator can 
be changed without changing the formalism of the operation itself. 

7) The same formalism can be used for pixel-based and other types of 
decomposition of the image (e.g. regions). 

[00133] In order to validate the method, we resolve the equation for 

steady state heat transfer in two applications: the linear diffusion and optical 
flow. These problems generate equations of the form of Equation 12 or its 
global version. We also use our methodology to resolve unsteady heat transfer 
with no source and we apply it to non-linear diffusion. 

Physical Principles (explanation of the physical foundations of the heat 
transfer problem) 

[00134] We distinguish two interesting particular cases for diffusion 

and optical flow problems: steady-state heat transfer and unsteady heat 
transfer with no source. 

[00135] Two important classes of laws are present conservation and 

constitutive laws. Conservation laws are independent of the properties of the 
material, whereas constitutive laws are specific to them. The physical 
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properties associated with a moving object are energy, work and heat In what 
follows, we describe each of these quantities. 

Energy Modeling 

[00136] Let us define some quantities for a continuous body 

occupying a volume V bordered by a surface S in a 3D space. We can say that 
such a body is composed of infinitely many particles (as many as desired), 
particles being the smallest elements [33]. Figure 5(a) illustrates such a body. 
At time t, a particle labeled X occupies a specific position: 

* = (*(*),»(*),*(*)) 



[00137] Each particle can move in space, so a velocity vector is 

associated with it at time t: 

v*(X,0 = ^ = v(x,t) 

[00138] Physical quantities can be associated with a particle labeled 

X (material description) or a position x tn space (spatial description). For 
example, v*(X,t) is the material description of the velocity of particle X and v(x,t) 
is the spatial description of the particle located at position x. For our purpose, 
we use spatial descriptions to derive the heat transfer equation. Vector n(x,t) is 
the outward direction of the surface at point x. 

[00139] The mass A m of a small amount of volume A V of a body is 

a measure of its inertia (tendency to resist motion). We use the term mass 
density, p to denote the following quantity: 
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[00140] p (x,t) is thus the mass density of the particle located at x at 

time t 

[00141] Two kinds of energy are associated with a moving object: 

kinetic and internal energy. Kinetic energy is a measure of the state of motion 
of a body: the faster the body moves, the greater its kinetic energy [28]. 
Because it is a measure of inertia, kinetic energy also takes the mass into 
account. For a particle located at x at time t, the kinetic energy is thus defined 
as 

K(x,t) = ip(x,t) (v(x,t) • v(x,t)) 

where . is the dot product. To obtain the kinetic energy for the entire body at 
time t, K(x,t) is integrated over the volume V: 

K(t) =\fff P(*> 0 < v (*>*) " v <*> *)) W 
where dV is an infinitesimal amount of the volume V. 

[00142] Internal energy is a measure of the state of temperature of a 

body. The hotter the body, the greater its internal energy. At time t, each 
particle has an internal energy density e (x,t) associated with it. The internal 
energy density is proportional to the temperature of the particle T(x,t) with a 
material-specific heat constant c; that is, s(x,t) = cT(x,t). For the entire body, 
the total internal energy is integrated over the volume V: 



E(t) = fff Pi***) e(x,*)dV 



(13) 
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[00143] The total energy for the body can now be defined as K(t) + 

E(t). 

Work Modeling 

[00144] Let us suppose that a body is submitted to an external force 

fe(x,t) (e.g. a traction force) and an internal density force b(x.t) (e.g. gravity). 
Figure 5(b) presents the action of external and internal forces. We define work 
as energy transferred to a body by means of a force acting on the body [28]. 
Work is negative when the energy is transferred from the body. Suppose that 
F(x) is an internal or external force that is constant over time, acting on a 
particle located at x during an amount of time t. This force will produce a 
displacement of the particle to position x*. This displacement is Ax = Xi - x. 
The work w(F, x) done by this force during this time is 

w(F,x) = F(x)«Ax 
and the instantaneous power P(x, t) is: 

P(x,t) = Urn ^Ei^ - F(x) • lim ^ = F(x) . ^ - F(x) • v(x,i) 

[00145] External forces act essentially on the surface of the body. 

The instantaneous work P e (t) done by the external forces on the entire body is 
thus the result of the integration of the external power over the surface: 

Pe(*) = y^^(x,*).v(x,t)dS 

where dS is an infinitesimal part of the surface of the body. 

[00146] Defining b(x, t) as the internal density force, the internal force 

is thus p (x, t)b(x, t). The rate of work over time done by the internal forces on 
the entire body is obtained by integrating internal work over the volume: 
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Pi(t) » JJJ p(x, t) (b(x, t) ■ v(x, *)) dV' 

The total work is thus P(t) ■ P e (t) + Pi(t) 
Heat Modeling 

[00147] Heat can be defined as energy transferred to a body owing to 

a difference in temperature. The heat flow density vector q(x, t) is a measure 
of the rate of heat conducted into the body per unit area per unit of time. We 
will see later how q(x, t) is defined. The external heat addition rate over time is 
the amount of heat coming from outside the body and entering by its surface. It 
is computed by projecting q(x, t) onto the inward normal vector (-n(x, t)) and 
integrating this projection over the surface: 

Qe(*)= ff q(x,t)-(-n(x,t))dS 

JJs (14) 

[00148] Now if a body has a rate of heat generation per unit of 

volume and time r(x, t), the internal rate of heat addition over time is computed 
by integrating r(x, t) over the volume: 

Qi(t) = JJJ P(x, t) r (x, t) dV 

[00149] Figure 5(c) shows q(x, t) and r(x f t). To simplify, we define 

the source a (x, t) as the rate of heat generated in a particle located at x per 
unit of volume and time: 



<r(x,*) =p(x,t)r(x 5 t) 



(15) 
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[00150] In many cases, this source Is known. However, it could also 

be a linear function of the temperature: <r(x, t) = a(x, t) + b(x, t)T(x, t) [38]. 
The total rate of heat addition over time is thus: 

Q(t)=Qe(t) + Qi(t) 

Energy Conservation Law 

[00151] A class of equations in continuum mechanics are those 

describing the conservation (equilibrium) principles. They express the 
conservation of certain physical quantities (mass, momentum, energy, etc.) 
over an entire body and as such they take the form of global equations over the 
whole body or a part of it [33], 

[00152] Conservation principles can be seen intuitively as follows: the 

change in the total amount of a physical quantity inside a body is equal to the 
amount of this quantity entering or leaving the body (through the boundary) and 
the amount generated or absorbed within the body. These laws are applicable 
for all continuous materials, moving and stationary, deformable and non-* 
deformable, and must always be satisfied. The global conservation equations 
can then be used to derive their local counterparts, called the associated field 
equations, which are valid at each point of the body including its borders. 

[00153] The first law of thermodynamics, which is relevant for the 

understanding of the heat transfer equation will now be discussed. This taw 
involves both kinetic and internal energies and states that the total variation of 
energy in a body (or a part of a body) is the result of the time rate of work and 
the rate of heat addition combined: 



%{E(t) + K(t)) = P{t)+Q(t). 



(16) 
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[00154] For heat transfer, we are only interested in the case of 

immobile bodies, that is v « (0,0,0), n(x,t) = n(x) and p(x,t) - /?(x). Equation 

16 thus becomes 

jijjj *) dV = J J -(q(x, t) • n(x,t)) dS + f ff^ a(x,t) rfV 

which now states that the thermal energy variation in a body is due to internal 
heat production added to the heat flowing into the body. Using the divergence 
theorem for Q e [33] and recalling that s(x,t) = cT(x,t), we obtain the thermal 
energy conservation law: 

where V . is the divergence operator. To simplify, let us define the temperature 

variation h(x,t) = — T(x,t). Equation 17 Is a conservation equation and is thus 

at 

valid over the entire body, a part or a point of this body. Consequently, the 
integral signs can be taken off: 

cp(x)ft(x,i) = — V-q(x,t) + cr(x,t). 

Thermal energy variation Rate of heat entering Rate of heat generation (18) 

[00155] This equation is said to be local, whereas Equation 17 is said 

to be global The thermal energy variation is called the unsteady term, the rate 
of heat entering is called the diffusion term and the rate of heat generation is 
called the source term. 

[00156] Based on Equation 18, two cases are be considered: the 

steady state case and the unsteady case with no source. The term steady 
simply means that there is no variation of the thermal energy of the system 
over time. That is, the left side of Equation 18 is null: 
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cp(x)/i(x,*) = 0 



(19) 



[00157] This implies that the heat diffusion compensates for the 

internal heat production: 

V-q(x,t) = a(x,t) (2Q) 
[001 58] In the unsteady case with no source we have 

<Kx,t) = o m 

which means that the time variation of thermal energy is explained by the heat 
diffusion alone: 

cp(x)/i(x, *) « - V ■ q(x, *) . 
Constitutive Principles 

[00159] In Equation 18, there are three unknown variables: p(x), 

h(x,t) and q(x,t). Let's look at the example of q(x,t). Suppose that we can 
measure the time variation of the thermal energy (left side of the equation) and 
also of the temperature T. We know that q(x,t) is related to the temperature, 
but since different materials usually have different diffusion properties, the 
missing equation q(x,t) = f(T,x,t) must depend on properties of the material we 
are studying, such as its homogeneity and type of diffusivity. Consequently, 
the system of equations contains more unknown variables than equations and 
the function f(T f x,t) must be added to the system formed by Equation 18. This 
is due to the difference in material properties. Different materials behave 
differently, but are subject to the same conservation laws. Constitutive 
equations such as f(T,x,t), which reflect the internal constitution of materials, 
allow us to complete the system of equations. 
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Decomposition into Basic Laws 

[00160] We have seen hereinabove that conservation equations are 

always valid regardless of the materials, whereas constitutive equations are 
dependent on their properties. When solving directly PDEs like Equation 18 in 
a discretized context with methods such as the finite differences approach, one 
makes global assumptions about the time and space behavior of the diffusion, 
energy variation and source terms without taking into account the nature of the 
basic laws underlying the problem. Some of these do not require any 
approximation since they come from conservation principles. Also, a more 
physically realistic solution can be obtained by choosing a proper 
approximation for each basic law arising from a constitutive principle. 
Consequently, we propose to decompose the terms of Equation 18 into basic 
laws. This equation can be broken down into one constitutive and two 
conservative laws for the steady state case. In the unsteady case, an additional 
constitutive and another conservative law must also be considered. Note that 
since the source term is often known, we do not try to decompose it. Recalling 
that the diffusion term a (x,t) is the rate of heat entering the particle located at x 
at time t, then 

a(x,t) = -V-q(x,t) 
is a first basic conservation law. 

[00161] The second conservation law concerns the thermal tension. 

We first define the thermal tension vector g(x,t) as the vector representing the 
direction and magnitude of the greatest temperature decrease at a fixed time t 
As g(x,t) is source-oriented (from hot to cold), we must put a minus sign before 
VT(x.t) which represents the direction and magnitude of the greatest 
temperature increase: 



g(x,*) = -VT(x,t) 



(24) 
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[00162] This equation is a second basic law. Since the thermal 

tension is the gradient of a scalar field, it is by definition a conservative field in 
space. We say that -T(x,t) is the potential field of g(x f t) [26, 41]. 

[00163] The third law is a constitutive law. We define the heat flow 

density q(x,t) as the quantity and the direction of the heat flowing into the 
particle located at point x at time t. It is represented by a vector and greatly 
depends on the behavior of the material. In the case of uniform, homogeneous 
materials, it has been proven experimentally by Fourier [20, 27] that q(x,t) is 
directly proportional to the difference of temperature relative to neighbors of this 
particle: 

q(x,*) = Ag(x,*) (2g) 

where X is a material-specific thermal conductivity constant. The value of X is 
known for many types of materials. Equation 25 is called the Fourier heat 
conduction law. For a non-homogeneous material, we consider that it has the 
behavior of a homogeneous material on an infinitesimal patch, but the 
conductivity changes with each patch; that is, 

q(x,*) = A(x,i)g(x,t) 

[00164] For the unsteady case, the fourth basic law is 

e(x,i) = cp(x)/i(x,t) 

where e (x,t) is the thermal energy variation (the unsteady term). This equation 
is a constitutive one because it involves p (x), which is material-dependent. 

[00165] Finally, the fifth basic law is related to the temperature 

variation and is a conservation law: 
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fc(x,*) = ^T(x,i) 

at (27) 

[00166] Considering only the temperature T x (t) of the particle located 

at x, we reduce the basic law to a 1 -dimensional equation and we can thus say 
that T x (t) is a conservative field in time-space. 

[00167] To summarize, we have three basic laws for the diffusion 

term of equation 18, that is: 

cv(x,t) = -Vq(x,i) 
q(x,<) - Ag(x,t) 
g(x,t) = -VT(x,t) 



[00168] We also have two additional basic laws for the unsteady 

term, that is: 

e(x,i) = cp(x)/i(x,t) 
h(x,t) = ^T(x,t) 



[00169] 



Combining all these elements, we obtain: 
cp(x, t)^-T{x, t) = V • (AVr(x, t)) + <r(x, t) 



(28) 



m 
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Discrete Representation of Images 

[00170] Some algebraic tools used to model images will now be 

recalled from the above-description. An image is composed of two distinctive 
parts: the image support (pixels) and some field quantity associated with each 
pixel. This quantity may be scalar (e.g. gray level), vectorial (e.g. color, 
multispectral, optical flow) or tensorial (e.g. Hessian). We model the image 
support in terms of cubical complexes, chains and boundaries. With these 
concepts, we are able to give a formal description of an image support of any 
dimension. For quantities, we introduce the concept of cochains, which are 
representations of fields over a cubical complex. For the use of these concepts 
in image processing, see [16]. 

[00171] As discussed hereinabove, an image is a complex of unit 

cubes usually called pixels. A pixel y c 9t" is a product 

7 = /iX/ a x,..x/ n 

where lj is either a singleton or an interval of unit length with integer endpoints. 
Thus lj is either the singleton {k} and is said to be a degenerate interval, or the 
closed interval [k, k+1] for some AeZ. The number q e{0,1 n} of non- 
degenerate intervals is by definition the dimension of y , which is called a q- 
pixel. Figure 6 illustrates three elementary pixels in 9t 2 . For g>\ 9 let 
J = {* 0 ,A,,..A^} be the ordered subset {1,2, ...,n} of indices for which 
I h =[a J9 bj] is non-degenerate. Let us define 

Akj<T = Ii X ... Ity-l X {dj} X I kj + l X ... X I n 

and 



B^O = hX ... hj-i X {bj} X x ... X I n 
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[00172] The A kj and the B kj are called the (q-1)-faces of a . One can 

define the (q-2)-faces, down to the 0-faces of a in the same way. The faces 
of y different from y itself are called its proper faces. 

[00173] By definition, a natural orientation of the cube is assumed for 

each pixel. Suppose that y denotes a particular positively oriented q-pixel. it is 
natural to denote the same pixel with opposite orientation by -y. Examples of 
orientations are given in Figure 6. A cubical complex in 9Tis a finite collection 
K of q-pixels such that every face of any pixel of the image support is also a 
pixel in K and the intersection of any two pixels of K is either empty or a face of 
each of them. For example, traditional 2D image models only considered pixels 
as 2D square elements. The definitions presented above allow us to consider 
2-pixels (square elements), 1 -pixels (line elements) and 0-pixels (point 
elements) simultaneously. 

[00174] In order to write the image support in algebraic form, we 

introduce the concept of chains. Any set of oriented q-pixels of a cubical 
complex can be written in algebraic form by attributing to them the coefficient 
0,1 or -1 , if they are not in the set or if they should or should not be taken with 
positive orientation, respectively. In order to represent weighted domains, we 
must allow arbitrary integer multiplicity for each q-pixel. 

[00175] Given a topological space X <zW in terms of a cubical 

complex, we get a free abelian group C q (X) generated by all the q-pixels. The 
elements of this group are called q-chains and they are formal linear 
combinations of q-pixels [16]. A formal expression for a q-chain Cq is 
c g-2L*Vi where A, eZ. 



[00176] The last step needed for the description of the image plane is 

the introduction of the concept of a boundary of a chain. Given a q-pixel y , we 
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define its boundary By as the (q-1)-chain corresponding to the alternating sum 
of its (q-1)-faces. The sum is taken according to the orientation of the (q-1)- 
faces with respect to the orientation of the q-pixeL We say that a (q-1)-face of 
7" is positively oriented relative to the orientation of yif its orientation is 
compatible with the orientation of y - By linearity, the extension of the definition 
of boundary to arbitrary q-chains is easy. For instance, in Figures 6(b) and 6(c), 
the boundary of the 1 -pixel a is X2 - Xi and the boundary of the 2-pixel A is 
we say that a and b are positively oriented with respect to the 
orientation of A but c and d are negatively oriented with respect to the 
orientation of A. Let us notice that the boundary of a 1 -pixel is always the 
difference between its boundary points. The boundary can be defined 
recursively. Given a (q-1 )-chain and a q-chain ^defined asy q = y^ x *[a 9 b] t the 
boundary of y q can be recursively written as 

37* = 97«-i x [a, b) + (-1)^(7^1 x {6} - y q _ x x {a}) (2g) 

[001771 In order to model the pixel quantity over the image plane, 

we look for an application F which associates a global quantity with all q-pixels 
yof a cubical complex. We denote this by <F>y>. This quantity may be any 
mathematical entity such as a scalar, a vector, etc. For two adjacent q-pixels 
y,and y 29 F must satisfy <F,Vi+V 2 >=4 <F 9 y l >+A 2 <F y y 2 >, which 
means that the sum of the quantity over each pixel is equal to the quantity over 
the two pixels. The resulting transformation F:C 9 (Jf)->9t is called a q- 
cochain and is used as a representation of a quantity over the cubical complex. 

[00178] Finally we need an operator which associates a global 

quantity with the (q+1)-pixels according to the global quantities given on their q- 
faces. Given a q-cochain F, we define an operator a, called the coboundary 
operator, which transforms F into a (q+1)-cochain 3F such that 



* • 
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<^,7>=< r,di> 

(30) 

for all (q+1>chains y. The coboundary is defined as the signed sum of the 
physical quantities associated with the q-faces of y. The sum is taken 
according to the relative orientation of the q-faces of the (q+1)-pixels of y with 
respect to the orientation of the pixels. Figure 7 presents an example of the 
coboundary operation for a 2-pixel. 

[00179] With this image model in hand, we will use the basic laws 

described hereinabove to rewrite the global heat transfer equation in algebraic 
terms [43]. 

Representation of the Heat Transfer Equation 

[00180] The process for representing the heat transfer equation in 

terms of algebraic topology can be summarized as follows. The image support 
is subdivided into cubical complexes. Basic laws are applied to pixels of various 
dimensions. These laws involve the computation of global quantities on pixels, 
expressed as cochains. Some of these laws link global quantities on pixels with 
global quantities on their boundaries and hence are expressed as 
coboundaries. The other laws are expressed as linear transformations 
between pairs of cochains. The topological formalism of cochain and 
coboundary is a generic one; that is, it does not offer computational rules. The 
cochains must be instantiated depending on the problem to be considered. 

[00181] The basic laws presented hereinabove will be reformulated in 

a topological way and then to give computational rules for cochains in the 
context of the heat transfer problem. Since we want to represent two kinds of 
global values over the spatio-temporal image, we use two complexes. The first 
complex is associated with global values corresponding to projections onto the 
tangential part of the domain (e.g. global thermal tension) while the second 
complex refers to values related to projections onto the normal part of the 
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domain (e.g. heat entering a particle). These two distinct orientations (see 
Figures 8(a) and 8(b)) give rise to two different complexes. 

Global Heat Transfer 

[00182] Let us assume that we have an image with n spatial 

dimensions and r pixels. Suppose also that we have a time interval [to, ti] which 
we can split into I equal sub-Intervals fc, V,], k e [0, 1-1]. Let us consider an n- 

complex representing the subdivided spatial support of the image K' . We can 
consider an (n+1)-complex representing the spatio-temporal support of the 
image: 

iC s = AC s/ x[^,4 + i],Vfce[0,Z-l] 

(00183] Now, let us consider y B% an (n+1)-pixel of K s , and the 

following cochain, defined as <e,y E >. Thus, we therefore need to define 
which value to use as cochain e in the heat transfer problem. Let us define s 
as the global energy variation on the (n+1)-pixel y E : 



< £,1e >= / e(x,t)dsy B 



(31) 



where dy £ is an infinitesimal part of the domain represented by y B . Now, 
using the global version of Equation 18, we obtain : 

/ c(x, t) dy B « / a (x, t) dy B + I a(x, t) dy E 

J 1B J tB JlE 



[00184] From this equation, we define two more cochains, 

representing first, the global diffusion: 
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(32) 

and second, the global source: 

< S, j E >= / cr(x, t) dy B 

[00185J Thus, we have the following relation between the three 

cochains: 

< £,7fi >=< T> t iB > + < S,j B > 

(33) 

[00186] The rules used for cochains sand D are then 

decomposed into basic laws. The rule for cochain S is not decomposed since 
we assume that its global value is known on y E . Let us finally mention that 
both steady and unsteady heat transfer problems can be considered using 
Equation 33 by setting respectively, cochains e and S, to zero. 

Global Temperature Variation 

[00187] Let us consider another n-complex, K" , representing the 

subdivided spatial domain of the image. We can define an (n+1)-complex 
representing the spatio-temporal image: 

K? = K?'x [t k) t k+l ),\/k € [0,/ - 1] 

[00188] Let us consider r „, a 1-pixel of K» defined as x^foAi], 

ie[1,rl. ke[0.Ml where xj is a 0-pixel of K"' . Let us also consider a 0-cochain T 
and a 1 -cochain H such that 
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< ^jh >=< 8T,jh >=< T,dy H > % 

(34) 

[00189] Figure 9 presents examples of T and H for K"of dimension 

3. 

[00190} Applying Equation 29, we find that the boundary of y„ is 

dY H - xix{tn + i} - Xix{tk}. According to the linearity of the cochain, the 
computational rule relating the global value associated to y„ with the values at 
its boundary X|x{t«} and Xjx{tK +1 } is 

< T,&y H >=< r,Xi x - ^ x {<*} >=< r,x { x {t k+i } > - < r,x, x {t k y > 

(35) 

[00191] This equation is general and applies to many problems. To 

define which values to use as 0-cochain and 1-cochain, let us take the global 
version of Equation 27 on y H and apply the fundamental calculus theorem: 

^ A(x, i) d 7 » = J j t T{x it t) dt = T(x f , t k+x ) - T(x i5 t k ) 

[00192] Looking at this equation, we see that it is similar to Equation 

35. Thus we define T = T(x, t). The location of the unknown temperatures we 
want to compute must correspond to the 0-pixels of K p . In order to fulfill 
Equation 34, we must therefore have 

<n,y H >= I h{x,t)dytt 

JyH (36) 

which is called the global temperature variation. These three equations are 
extended by linearity to a 1-chain of defined as yxfo, where yis an 
arbitrary 0-chain of K p . 
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Global Energy Variation 

[00193] We want to link cochains H and e, representing the giobal 

temperature variation and the global energy variation, respectively. For this 
purpose, we need to represent Equation 26. The two cochains are not from the 
same cubical complex (H is from tf'and e is from K s ), and moreover, 
Equation 26 is material-dependent; therefore they cannot be linked exactly. 
However, we can express this link as a linear transformation 

H— z-+£ 

[00194] Recalling Equation 31 , we have 



£,7£?>= / €(x,t)d 7£ = / cp(x)fc(x,r) dy B 



(37) 



[00195] Unfortunately, we do not know the value of p (x) or h(x, t) at 

all points of the volume. Consequently we must approximate these two fields 
over the volume. We denote the approximations by p(x) and £(x. t). For one 
1-pixel y H , defined as xixflk.tk+i], the approximation is performed piecewise 
such that h(x, t) must fulfill Equation 36. 



Mx*,t)s<W,7£r > 

(38) 



[00196] Equation 37 thus becomes 

< €, y E >= / c/3(x)ft(x, t) d lB = /e(c, H) = T 

J " B (39) 

where dV is an infinitesimal part of y E which depends on the choice of the 
approximation functions p{x) and h(x, t) and the position of K * with respect to 
K". 
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Glofoal Diffusion 

[00197] Let us consider an n-cochain Q and an (n+1>cochain D 

defined by the coboundary 

< 2>, 7b >=< <5Q, 7b >=:< Q, dy E > ( 

[00198] Figure 1 0 presents examples of Q and D for K' of dimension 

3. Let us assume that the n-faces y a of y E are positively oriented relative to 
y B . According to the linearity of the cochain, the computational rule is 

<Z>,7e>= J2 <Q,7o ( > 

(41) 

[00199] Again, this equation is general; hence we need to find a 

global value for the (n+1>cochain D, which can be computed by summing the 
global values at the boundary of y E . According to Equation 32, we have 

<2?, 7B >= f a(x,t)dy B 
[00200] Applying the divergence theorem to this equation, we obtain 

where n(x, t) is the normal vector to an infinitesimal part of the domain 
represented by y e . This last equation is in the form of a coboundary (Equation 
41), from which we define 

< Q, 70i >= / -q(x, t) • n(x, t) dy Qi 

(42) 



CA 02397389 2002-08-09 



-48- 

[00201] Again, the previous definitions can be extended by linearity to 

arbitrary (n+1)-chains of K' . We recall that there is absolutely no 
approximation in these equations. 

Global Thermal Tension 

[00202] Let us consider a 1 -pixel y G of defined as yxt k , 

k e[0,/-i], where y 'S a 1 -pixel of K" whose boundary is defined as dy = x, - 
xi, i,je[l,r]. Let us also considera 0-cochain T and a 1-cochain G defined by 
the coboundary 



< 0, 7c >=< <JT, 7o >=< T, djo > {43) 

[00203] Figure 1 1 presents examples of T and G for "one 3-pixel of 

K' . The boundary of y 0 is dy a * xj x{tk} - xi xft}. According to the linearity of 
the cochain, the computational rule relating the global value associated with 
y a to the values at xi x{tk} and xj x{tk} is 



< T,dy a >=< T,xj x {t k } - x {t k } >=< r,x, x {<*} > - < r,Xi x {t k } > 

(44) 

[00204] To define which values to use as cochains G and T let us 

take the global form of Equation 24 on y G 

f S(M) • rfTG = f g(x, t k ) - dy = p -Vr(x, i fc ) • rfy 

J no J Xi J* (45) 

where dy is an infinitesimal part of y. Since g(x, t) is a spatial conservative 
field, we can apply the line integral theorem [26, 41] saying that for a 
conservative field F(x) = Vf(x) and two points A and B, in an open connected 
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region containing F(x), the integral of the tangential part of F(x) along the curve 
R joining A and B is independent of the path (Figure 12): 

f P(x) • dR x = f* P(x) • dRs = f B P(x) • dR 3 = /(£) - f{A) 
j a Ja Ja 

[00205] From this theorem, we can rewrite Equation 45 as 

/ ' g(x, h) -dy = (-T( Xj , < fc )) - (-T( Xi , = r(xi, t fc ) - T(x is t fc ) 
Jx< (46) 

[00206] Looking at Equation 46, we see that it is similar to 

Equation 44. Thus we define T » T(x, t). Consequently, the location of the 
unknown temperatures we want to compute must correspond to the 0-pixels of 
K" which is coherent with the conclusions hereinabove. In order to fulfill 
Equation 43, we must have 

< G, 1g >= I ~g(x, t) • dy G 

J ™ (47) 

[00207] The previous definitions are extended by linearity to 1 -chains 

of K" defined as y x^, where y is an arbitrary 1 -chain of K"' . 

Heat Flow Density 

[00208] The coboundaries <Q > dy E > (Equation 40) and <T,dy c > 

(Equation 43) provide exact global versions of Equation 23 on K' and Equation 
24 on K" , respectively. In order to complete the diffusion term, we need to 
represent Equation 25, which links local values g(x,t) and q(x,t). Equation 25 is 
a constitutive equation and cannot be represented by a topological equation. 
However, we can find a relation transforming cochain G into cochain Q: 



<£,7<7> A <Q,7Q> 



CA 02397389 2002-08-09 



-50- 

as a global counterpart for Equation 25. To find this transformation, we recall 
Equation 42: 

< Q» 7<?, >= / -q(x, t) ■ n(x, t) dry Qi = / - Ag(x, i) • n(x, t) rf 7<?( 

relating cochain Q to field g(x,t). Unfortunately, field g(x,t) is not known, so we 
have to approximate it with a field g(x,t). Let us consider y a , an n-pixel of 

defined as r„=r, *e[u,/-l] where y x is an n-pixel of K"' . This 
approximation is performed piecewise such that for each 1-face y c of y„ , g (x,t) 
satisfies 

Jn,a (48) 

where dR is an infinitesimal part of the domain represented by y G . Equation 25 
is then applied to obtain g (x.t): 

q(x,t) = Ag(x,*) 
at all points of the domain. Equation 42 becomes 

< QtfQt >= / -q(x, t) • n(x, t) d lQi = j g {\ Q) 

J " Q * (49) 

[00209] The transformation we are looking for is thus 

A = / fl (A,S) 

which depends on the choice of an approximation function g(x,t) and the 
position of K' with respect to K" . 



CA 02397389 2002-08-09 



-51- 

Boundary Conditions 

[00210] The decomposition process we have presented herein is 

carried out with the assumption that all the needed quantities surrounding a 
pixel are known. For instance, in the steady state heat transfer problem, for a 
particular (n+1)-pixel, the cochain S must be known for all other surrounding 
(n+1)-pixels, that is, there are as many equations as variables. 

[0021 1] Unfortunately, this assumption is not verified at the borders of 

the image. Thus, as in solving the PDE, certain boundary conditions must be 
imposed to specify the gray-level conditions at the boundary of the image. For 
instance, these conditions may prescribe the values of either cochain T 
(Dirichlet boundary conditions) or cochain Q (Neumann boundary conditions). 

Summary of the Algorithm 

[00212] We will now summarize the algorithm used to find an 

expression of the temperatures at time t* +1 as a function of the temperatures at 
time tk. The input data for this algorithm are the cochain S and the Dirichlet 
boundary conditions. That is, T is known for all pixels on the boundary of the 
image, which includes the values at time to. 

1 . Choose the positions for K p> and K' . 

2. Compute e as a function of H: 

(a) Choose the approximation functions h (x.t) and p (x). 

(b) Apply Equation 38 to find h (xA, t,) as a function of H. 

(c) Apply Equation 39 to find the transformation r, expressing s as a 
function of H. 



3. Apply r and Equation 34 to find s as a function of T. 
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4. Compute Q as a function of G: 

(a) Choose the approximation function g (x,t). 

(b) Apply Equation 48 to find g (x,t) as a function of G. 

(c) Apply Equation 49 to find the transformation A , expressing Q as a 
function of G. 

5. Apply Equation 40, A and Equation 43 to find D as a function of T. 

6. Apply Equation 33 to obtain an equation of the temperatures at time tfe + i as 
a function of the temperatures at time tk 

[00213] Figure 13 presents an overview of the computational scheme. 

Applications 

Linear isotropic Diffusion 

[00214] One of the most direct applications of the heat transfer 

equation is the isotropic diffusion of gray-level intensities; that is, smoothing. 
For a 2D image l(x), with x * (x,y), the resolution of the PDE 



dt 



J(x,*) = V 2 I(x,t) 



(50) 



is equivalent to the convolution 



/(x,*) = (J*&)(x) 



where 
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is a Gaussian with variance a 2 ~2t\25\. One can see t as the scale of the 
smoothing operation. Let us assume that the Laplacian image at scale t 

L(x,*)=V 2 J(x,t) (51) 

is known. One can consider this equation as a steady state heat transfer 
problem with T(x,t) = Kx.t), <x(x.t) = -L(x,t) and A = 1. 

[00215] We want to solve Equation 51 for local l(x,t) located at the 

center of each image pixel. Employing the process presented hereinabove, we 
first position the two cubical complexes representing two subdivisions of the 

image plane. As stated hereinabove, the primary complex K"' is defined with 0- 
pixels corresponding to pixel centers. For the sake of simplicity, 

K' corresponds to the image pixels; that is, the secondary 2-pixels y, are 
rectangular and symmetrically staggered relative to the 1 -pixels of K" and the 
1 -pixels y q of K' intersect orthogonally in the centers of the primary 1 -pixels. 
Since there is no variation in steady-state heat transfer over time, we drop the 

time parameter. This means that K P = K>, = and the time integral in 
cochain computation are dropped. We saw that the approximation function f a 
depends on the position of K' with respect to K" . 

[00216] Figure 14 shows the two complexes for a 5 x 5 image. 

Positioning the 1 -faces of K* such as each passes through the center point 
between two 0-pixels of K" allows us to compute a polynomial function of order 
1 with the same accuracy as that obtained using one of order 2 [37, 34]. 

[00217] We need a global value for the 2-cochain <S,y E >. If we 

assume that a pixel value represents the global value of intensity, we can 
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directly set <S,y E > = -L(x). This assumption is reasonable if we look at image 
acquisition as a process which accumulates the total number of photons within 
a global area corresponding to the pixel [221. 

[00218] We must choose an approximation function g(x). For 

simplicity, we assume that g (x) arises from a bilinear approximation; that is: 

g(x) = (a+by) ' (c + dx) • j 

[00219] Given a 2-pixel y p of K" , g (x) must satisfy Equation 48 

for each 1-face of y p . As an example, let us find the coefficients a, b. c and d 
for such a pixel defined as in Figure 15. We have 

/•a 

Gi — I -g(x, 0) • idx 
Jo 

02 = f -g(A,y).jrfy 
Jo 

£?3 = f -g(x, A) • idx 
Jo 

04 = f -g(0,y) Jdy 
Jo 

from which we obtain 

•« - 4 [(* - ■*+ (ft- • J] . « * w 

£(x) is thus a piecewise function of G, but as G is computed from T, we can 
also express g(x) as a function of T. For each primary 2-pixel, we apply 
Equation 25 to obtain g(x)~ g (x). 
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[00220] The next step is to compute <Q,y Q > for K' from 

Equation 49. Each secondary 2-pixel y M intersects with four primary 2-pixels, 
Y pa . 7>. Ypc and Yn • There are four segments in the approximation function 
q (x) corresponding to the four primary 2-pixels; that is, q a (x), q b (x), £(x) and 
?„(x). Figure 16 illustrates y E . We find cochain <Q,y Q > corresponding to the 
four 1 -faces of y E : 

' -q«( A/2, #) -Tdy+ -q 6 (A/2, j,) • «fy 

0 •/ —A/2 

3(7t & , (?5 

T~ Y T 

^ -q 6 (x, -A/2) • {-j)dx + y -q c (a;, -A/2) • (-J)da 



4 8 8 



03 = ^ -qa(a:, A/2). Jdx + y -q d (ar, A/2) • Jdx 



3^ , & _,_ 0n 



24 = 7„ -^(-A/2 l2 /)-(-i)dj,+ / -q c (-A/2,y).(-?)dy 

* 0 •» —A/2 



~ 4 8 ~~8~ 

(53) 

[00221] Using Equation 41 , we have 
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<V,7B>= Q1 + Q2 + Q3 + Q4 tKA . 

(54) 

[00222] Substituting Equation 43 in Equation 53, Equation 53 in 

Equation 54 and Equation 54 in Equation 33, we can now express <S,r B > as 
a function of T . As an example, we present <S,r E > for a 2-pixel y E , defined 
as in Figure 16: 

< 5, 7a >= -373,o + 1 [75,! + T U o + To,-* + 71 1>0 ] + \ [T-i., + T ul + T^i + T-t.-t] 

(55) 

[00223] For each non-border pixel (represented by a secondary 2- 

pixel), we get an equation in the form of Equation 55. For the border pixels, we 
set T = l(x). Solving this system, we obtain the smoothed image l(x,t) « T. 

Optical Flow 

(00224] An indirect application of the heat transfer equation is the 

computation of optical flow for a 2D image sequence l(x,t), using the Horn and 
Schunk [29] algorithm. It can be shown that the velocity vector u(x,t) = (u(x,t), 
v(x,t» satisfies the following constraint arising from variational calculus (for 
greater legibility, we have dropped the (x.t)) : 

I 2 u + I x J y v = a 2 V 2 u-I x I t 
Ixlyu + ftv = a 2 V 2 v - I v I t 

v v (56) 

where a is a weighting factor and l x , l y and l t are the first derivatives of l(x,t) in 
x, y and t, respectively. Let us rewrite Equation 56 in the following vectorial 
form: 



V/(V/ • u) = or 2 V 2 u - J t VJ 
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where we define V 2 u = (V 2 u, V 2 v). Reorganizing the terms of Equation 56, we 
get the following equation: 

a 2 V 2 u = VJ(VI • u) + / t V/. 

(Of) 

[00225J Taking <r(x, t) = -VI(V| . U ) - l t VI as a heat source, 

Equation 57 can be seen as a steady state heat transfer equation in which the 
cochain T corresponds to u(x, t) and X =a 2 . It can thus be decomposed using 
the method descrbed hereinabove. The cochain T is u(x, t), and we get the 
relation: 

[00226] For the same reasons as in the linear diffusion problem, 

special considerations are needed at the borders of the image. We assume 
zero velocity at the borders of the image and solve the system to get the 
velocity field for each point of the image. 

Nonlinear Diffusion 

[00227] Linear isotropic diffusion reduces noise but also bturs edges. 

As the scale increases, edges tend to be harder to identify [43]. One possible 
way of reducing this effect might be to consider the heat conduction coefficient 
X as a field function dependent on the magnitude of the edges; that is: 

Jf/(x,0 - V • (ff(|V/(x,i)| 2 )V/(x,*)) 

m (58) 

which corresponds to Equation 28 with A(x,t) - g(|VI(x f t)I 2 ) f T(x,t) = l(x,t), 
p(x) = 1, c 9 1 and <r(x,t) = 0 (i.e., unsteady transfer with no source). The 
conduction function g(s) must display the following behavior, in constant 
regions, there should be linear isotropic diffusion (Equation 50), that is 
g(fVl(x,t)| 2 ) = 1 for |VI(x,t)| 2 « 0, and almost no diffusion when the magnitude 
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of the edge is great; that is, g(|VI(x,t)| 2 ) = 0 for |VI(x,t)| 2 -»•«>. Perona and 
Malik [38] proposed the following functions: 

and 

g(s) = e **, (& > 0) 

[00228] The parameter k in these functions is difficult to set because 

it controls the threshold of diffusion but also the steepness of the function [35J. 
We prefer to use the function 

g{s) = |tanh (<y(k - s) + 1) 

where k and r control the threshold and the steepness, respectively. We then 
solve Equation 58 for a particular t (the scale) with initial conditions 

J(x,0) = J(x) 

where l(x) is the original image. 

[00229] Let us assume that we have I time steps Ar = t II . First, we 

use the same cubical complexes K" and K' as hereinabove and we define 

K' = /C a 'x t k = kAt, V* 6 [0,1-1] 

[00230] Secondly, we make an assumption about the spatial behavior 

of h(x, t); that is, we choose the approximation function h (x, t). For a 3-pixel y E 



CA 02397389 2002-08-09 



-59- 

as defined in Figure 17, let us assume that H is the mean value over 
[-A/2,A/2] x [-A/2, A/2], We thus have s = H ; that Is, using Equation 34, we 
can also write 

[00231] For the sake of simplicity, we use the same spatial bilinear 

approximation function g (x, t) as hereinabove. We have to approximate the 
behavior over a time step. Some common assumptions about time variation 
may be generalized by proposing [37]: 

/ A(t) dt = (wA(t k+1 ) 4- (1 - w)A(t k )) At, 0<w<l 

where A(t) is some quantity and w is a weighting factor [37]. Some values of w 
lead to well-known schemes: 1) w leads to the explicit scheme; that is, the 
value at % prevails for the entire time interval except at time W 2) w = 1 leads 
to the fully implicit scheme; that is, the value changes at time tk from A(tk) to 
A(tn +1 ) and stays there throughout the whole time interval. 3) w = 0.5 leads to 
the semi-implicit or Crank-Nicolson scheme; that is, there is a linear variation of 
A(t) . We propose to use the implicit scheme because for large values of At , it 
best emulates long term time behavior for heat [37]. That is for a 3-pixel, we 
have for w = 1 

tw - -h [(« + ■ (« + •/] x . . to, a.) 

[00232] In order to obtain the local function q (x. t) we apply Equation 

25: 



q(x,t) = A(x,*)g(x,t) 
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where A(x, t) = g(|VI(x,t)| 2 ). As Vl(x,t) is a spatially sampled Image where 
samples are located at the 0-pixels of K" , we must approximate the local 
values of A (x, t). For the sake of simplicity, we once again use a bilinear 
approximation; that is: 



A(x) = a + bx + cy + dxy 
For a 2-pixel of K" , as illustrated in Figure 18, we obtain 



[00233] 



A(x, t) = A 0 ,o + ((A t , 0 - A 0)0 )a; + (A 0 ,i - A 0f0 )y) + i(A 0 ,o - A 1>0 - A 0 ,i + \i±)xy 



A 2 ' 



(60) 



[00234] Using these assumptions, we follow the same steps as in 

hereinabove to find Q as a function of G. For instance, this function for one n- 
pixel as defined in Figure 16(c) is 

Qt = (CL) l G 
where C, L and G are matrices defined as 



-[ 



1/24 7/24 1/24 1/24 7/24 1/241 

0 1/24 1/48 0 1/24 1/4B , 
1/48 1/24 0 1/48 1/24 0 J 



L = 



Ao.-x" 
*o,o 
Ao.i 

*i,o 



and G 



■a 



[00235] Using Equations 40 and 43, we can express D as a function 

of T . For one 3-pixel, y B of K' as defined in Figure 19. we have 



<Z>,7i?>= (CaLa/TA* 
where C, L and T are matrices defined as 



(61) 
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1/24 1/16 0 

1/48 1/2 1/48 

0 1/16 1/24 

1/48 0 0 

-1/12 -3/8 -I/" 

0 o i/ 4 » 

Q O 0 

o o o 

L o u o 



1/16 1/12 0 

0 5/24 0 

0 1/12 1/16 

1/4 5/24 0 

-3/8 -7/6 -3/8 

0 5/24 1/4 

1/16 1/12 0 

0 5/24 0 

0 1/12 1/16 



0 
0 
0 
1/48 



0 
0 
0 
0 



0 o 

1/24 1/16 

1/48 1/4 

0 1/16 



0 
0 
0 
-1/1 
1/48 
0 



* 






T.-» 










A- 1,0 




T fP 

'0.0 

% 




Ao.o 


and Ta 


Xi.o 
A-M 
Ao.i 
. Ai.i 


1 



[00236] We apply Equation 33 with <S,y E > = 0 and we obtain, for 

9 

each 3-pixel of K' , 

7^ 0 -(C A L A ) t TAf = 7? 0 

which defines the system of linear equations. The initial conditions are T° = 
l(x). 

A Different Hypothesis for Heat Conduction 

[00237] In the preceding discussion, we approximated A(x) with a 

bilinear function, essentially for the sake of simplicity. Nevertheless, it could be 
preferable to use another assumption. Actually, this simple approach does not 
accurately handle abrupt changes in conductivity. For example, let us consider 
two 2-pixels of K', as shown in Figure 20. To compute Q, we need to 
approximate A(x) at the borders of the pixels based on the values at their 
centers. Using bilinear approximation, the value of X (x) on the line linking two 
points is declared be the arithmetic mean of the values at these points. For 
instance, given ^ -» 0 and \ 0 - 1 . the conduction at the border is about 0.5. 
This means that the zero conductivity at one pixel is partly cancelled out by the 
fact that on the pixel beside it. there is a high conductivity coefficient In non- 
linear gray level diffusion, we are confronted with precisely this kind of abrupt 
change. For example, at step edge pixels, the conduction may needs to be 
very low, whereas immediately adjacent, it may needs to be almost one. 
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[00238] A better assumption would thus be to consider X(x) as 

constant over one single 2-pixel of K s [37]. Therefore, on the 1-face common 
to two pixels as in Figure 20, we have 



\Ao,o Mfi/ Ao,o + Ai 



fl_ 
Ai,o 



[00239] It can easily be seen that when \ 0 -»0, then X-»o and 

when \ 0 «\ 0 , then X-+2\ 0 . This means that in both situations, the low 
conductivity would prevail at the boundary common to the two pixels {37J. With 
this assumption, the matrices C x and L x are modified as follows: 



-1/4 


1/4 


0 


o ■ 


3/2 


-1/4 


-1/4 


0 


1/4 


0 


1/4 


0 


-1/4 


3/2 


0 


-1/4 


-3/2 


-3/2 


-3/2 


-3/2 


-1/4 


0 


3/2 


-1/4 


0 


1/4 


0 


1/4 


0 


-1/4 


-1/4 


3/2 


. 0 


0 


1/4 


1/4 J 



, and h\ =5 Ao,o 



Ao.-i/(Ao,_i + Ao.o) 
A_i,o/(A_i 9 o + Ao.o) 

Ai f o/(Ai f o + Ao t o) 
, Aoa/(Ao f i + A 0 ,o) 



Experimental Results 

[00240] The proposed approach was tested on real and synthetic 

images in the context of linear isotropic diffusion, optical flow and non-linear 
diffusion. The results were compared with another method in each case. 

[00241] For linear diffusion, Figure 21(a) presents our physics-based 

method (PB) at three different scales and Figure 21(b) shows the result by 
convolution for the same scales. In the absence of a quantitative evaluation, 
we can say that subjectively the results seem to be similar which fulfills our 
objective of validating our approach. 
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P0242] 



For optical flow, Figure 22 shows the first frames of three 



sequences: the rotating sphere, Hamburg taxi and tree sequences. The results 
are compared with those being obtained using a finite-difference 
implementation of the Horn and Schunck algorithm (FD) [18, 19]. In these three 
examples and for both the PB and FD methods, the image derivatives are 
computed by convolution with the appropriate Gaussian derivatives. Both 
temporal and spatial scales are set to 1, as is the weighting factor a . Figure 23 
shows the flow pattern computed for the sphere sequence. Figures 24 and 25 
present the flow patterns for the taxi and tree sequences, respectively. For the 
rotating sphere and the taxi sequences, we obtain similar results with both 
methods. For the tree sequence, we also obtain similar results even if the 
extreme values seem to be smaller with the PB method than with the FD 
method. This fact is more apparent in Figure 27(a) and 27(b) where we show 
respectively the results for the PB and FD methods for the tree sequence in 
which white noise (standard deviation of 10) has been added (see Figure 26). 
Another advantage of our method is that it avoids iterations since we apply the 
algorithm only once. 

100243] For nonlinear diffusion, Figure 28 compares the PB with 

constant hypothesis on X and FD [38, 17] methods for a small window of the 
peppers image with a = 5. Figure 28(b) presents the original section with white 
noise added (standard deviation of 10). Figures 28(b) and 28(c) show 
respectively the results for the PB and FD methods. One can notice that some 
details are better conserved in Figure 28(c) than in Figure 28(b). This fact is 
enlightened in Figure 29 which shows a profile of the diagonal line starting at 
the upper right corner and finishing at the lower left corner. 

100244] Figure 30 presents the results for the peppers image with a 

= 1 .0, 3.0 and 5.0. The results for PB seem a little sharper than the FD results. 



[0024SJ Figure 32 shows the results for the Lena image (Figure 31(a)) 

with an added white noise of standard deviation 10 (Figure 31(b)) at two 
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different scales, a = 4.0 and 8.0. Again, the PB method seems to give sharper 
results at both scales. 

(00246] Figure 33 presents details of the Lena results at cr = 8.0. 

Figure (b) seems smoother in constant zones but some details are lost. For 
example, compare the eyes, the trim on the hat and the right of the face. 

[00247] An alternative approach to the PDE-based resolution of the 

diffusion problem was described. The proposed approach differs in two 
significant ways from with the classical PDE resolution scheme: 1) the image is 
considered as a cubical complex for which algebraic structures such as chains, 
cochains, boundaries and coboundaries are defined; and 2) the diffusion 
problem is decomposed into conservative and constitutive basic laws, each of 
which is represented by cochains and coboundaries. 

|00248| The conservative basic laws are represented without 

approximation while some approximations are required for the constitutive 
laws. This means that unlike traditional PDE resolution, for which many 
approximations must be made, all approximations are known since they are 
only needed in the representation of the constitutive equations. Coboundaries 
are computed using fundamental theorems of calculus such as the Green, 
Stokes and Line Integral theorems. Unlike iterative numerical analysis 
algorithms that do not allow the explanation of intermediate results, the use of 
basic laws allows the physical explanation of all steps and intermediate results 
of the algorithm. Moreover, since there is no iteration in the resolution process, 
there is no problem about the convergence of the numerical analysis scheme. 
Furthermore, the use of cubical complexes provides algorithms that can 
operate in any dimension. It has the significant advantage of avoiding the 
potentially difficult task of extending the algorithm to higher dimensions. 
Cochains and coboundaries allow the use of both global and local quantities. 
Integrals or discrete summations over fields are used to compute global 
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quantities. This allows the reduction of noise by performing a smoothing 
operation, as opposed to differentiation, which enhances high frequencies. 

[00249] In computer vision and image processing, several problems 

can be modeled as diffusion problems. The proposed approach has been 
validated on smoothing by linear and nonlinear diffusion and on the 
computation of optical flow. The results obtained confirm the effectiveness of 
this approach. 

PRACTICAL EXAMPLE #2: 

A PHYSICS-BASED MODEL FOR ACTIVE CONTOURS 

[00250] A new active contours model is presented. It is based upon a 

decomposition or the linear elasticity problem into basic physical laws. As 
opposite with other physics-based active contours model which solve the partial 
differential equation arising from the physical laws by some purely numerical 
techniques, we use exact global values and make approximations only when 
they are needed. Moreover, these approximations can be made wisely 
assuming some knowledge about the problem and the domain. The 
deformations computed with our approach have a physical interpretation. In 
addition, the deformed curves have some interesting physical properties such 
as the ability to recover their original shape when the external forces are 
removed. The physical laws are encoded using the computational algebraic 
topology based image model described herein. The resultant numerical 
scheme is then straightforward. The image model allows our algorithm to 
perform with either 2D or 3D problems. 

Introduction 

[002S1] These last years, active contours and active surfaces have 

been widely studied since the introduction of active contours by Kass et al [59]. 
They have been used in image segmentation [62], tracking [68], automatic 
correction and updating of road databases [46], etc. 
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£002521 To solve these problems, many different approaches have 

been proposed (see [57, 63]) in particular physical models derived from 
equations of continuum mechanics. Mass-springs models are physical models 
which use a discrete representation of the objects. Objects are modeled as a 
lattice with point masses linked together by springs [57]. information is thus 
only available at a finite number of points [63). These methods offer only a 
rough approximation of the phenomena happening in a body [57]. Moreover, 
the determination of springs constants reflecting the material properties may be 
a very fastidious work. However, they offer real-time performances and allows 
for parallel computations. 

£00203] Other physical models based upon the minimization of an 

energy functional which takes into account an internal regularizing force and an 
external force applied on the data are also often used. Some of them consider 
the deformable bodies as continuous objects by approximating their continuous 
behavior with methods such as the finite element method (FEM). FEM are 
closer to the physics than mass-springs models but their computational 
requirements make them difficult to be applied in real-time systems without 
preprocessing steps [57]. Finite difference methods (FDM) are also used to 
discretize the objects. They usually offer better performance than FEM but they 
require the computation of fourth order derivatives which make them sensitive 
to noise [63]. 

K002S43) For a given curve S, the application of the FEM and FDM 

methods leads to a discrete stationary system of equations of the form 

KS = t(S) 

where K is a matrix which encodes the regularizing constraints on S and 1F(S) 
represents the data potential. However, some problems such as animation in 
graphics applications require to take into account a dynamic evolution of the 
curve [57]. In these case, inertial body forces and damping forces may also by 
considered by controlling the deformations through a Newtonian law of motion 
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at 2 at (62) 
or by a Lagrangian evolution 

™+KS = f(S) 

at (63) 
where M and D are respectively matrices which represent the mass model and 
the background damping. Equations 62 and 63 are solved using various 
numerical schemes [70, 64] assuming an initial curve S 0 close to the solution 
which evolves until the inertial terms go to zero. 

[00255] Over the last years, a lot of different methods have been 

introduced to compute the matrices M, D and K but as pointed out by 
Montagnat et al [63], these methods have a major drawback: the corresponding 
system deformations do not have any physical interpretations. 

[00256] A new model which includes a physical interpretation of the 

deformations is described hereinbelow. The model is similar to a mass-springs 
model but it includes both the efficiency of the mass-springs models and the 
accuracy of the physical modeling of the FEM by providing a systematic 
method for specifying springs constants reflecting the properties of the 
materials. 

[00257] To achieve it, we propose to use directly the basic laws of 

physics which lead to the partial differential equations 62 and 63. These 
equations are indeed obtained by considering and mixing together some basic 
laws of physics into a global conservation law and considering its local 
counterpart [72]. This approach is not always well suited for problems such as 
computer vision in which the continuous domain must be subdivided into many 
sub-domains for which there are often only one information available. The use 
of this information as a global value over each sub-domain allow to directly use 
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the global conservation law which can lead to an algorithm less sensitive to 
noise. 

[00258] To encode these global values over points, surfaces, 

volumes, etc arising from some physical laws, we use the computational 
algebraic topology based image model described hereinabove. 

[00259] Our approach has several advantages. 1) Since the linear 

elasticity problem is well-known in continuum mechanics, our modeling can be 
made wisely in order to provide some good physical interpretation of the whole 
deformation process and of its intermediate steps. This allows an easier 
determination of the parameters used in the process since they have a physical 
meaning; 2) The determination of the springs constants in order to reflect the 
material properties is straightforward; 3) The objects in the image (e.g. curves, 
surfaces) are modeled as entities having their own physical properties such as 
elasticity and rigidity. They have the property of recovering their original state 
when the forces applied on them are removed; 4) Both smooth results and 
results having high curvature points can be obtained; 5) The complexity of the 
algorithm is minimal and allows for real-time simulation without any extra 
preprocessing steps [51]; 6) The image model allows our algorithm to perform 
with either 2D or 3D problems. 

Physical Modeling 

[00260] One of the objectives of this aspect of the present invention is 

to model the objects in an image (e.g. curves, surfaces, etc) as entities having 
their own physical properties such as elasticity and rigidity. As a consequence, 
these objects need to satisfy the laws and principles of the continuum 
mechanics. For instance, a body subjected to forces must move or deform 
according to the universal laws of physics. 

[00261] These principles and laws to which all bodies must obey will 

now be presented. We first introduce the concepts of stress and strain which 
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are essential in the statement of the governing equations for deformable 
bodies. We then present the physical laws related to the linear elasticity 
problem. 

[00262] The elasticity theory has been widely studied by engineers 

and scientists and is the main subject of many books. We only present the 
concepts of this theory which are relevant to our application. The concepts 
presented here are well-known and may be found in many continuum 
mechanics books such as [65, 50]. 

Forces, Stresses and Strains. 

[00263] A material body in a 3-D space Is always subjected to forces. 

These forces may come from an external agent (external forces) or issue from 
the object itself (internal forces). When the external forces are greater than the 
internal forces, then the body can undergo deformations (strains) or be 
accelerated. This deformation can induce internal forces (stresses) if the 
material is elastic. The concepts of force, stress and strain and the relation 
between strain and stress is exposed hereinbelow. 

Forces acting on a body 

[00264] Two basic types of forces act on a body. First, there are the 

interatomic forces which hold the body's particles together at some 
configuration. These forces, called internal forces, can either attempt to 
separate or bring the particles closer according to the fact that the body 
undergoes a contraction or an extension [65]. They act in response to a force 
applied by some other agent Assuming the equilibrium of the body and using 
Newton's law of reaction, they must be equal in magnitude to the forces applied 
on the body but in opposite direction [66]. 

[00265] On the other hand, there are the forces applied by an 

external agent called external forces. Two types of these forces are generally 
applied on a body. First, there are forces such as gravity and inertia called the 
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body forces, which act on ait volume elements. These forces, noted b t (forces 
per unit of mass in a direction Xj), are distributed in every part of the body. 
Secondly, there are forces which act on and are distributed over a surface 
element such as the contact forces between solid elements [49]. These forces, 
noted f, (forces per unit of area in a direction xj, are called the surface forces. 
The surface element may be inside the body or a part of a bounding surface 
[60]. A body of arbitrary size, shape and material subjected to surface and body 
forces is shown in Figure 34. 

[00266] External forces applied on a body must be transmitted to it. A 

rigid body can then undergo either a spatial shift, a rotation of both of them. In 
the case of a non-rigid body, it can also go through a deformation or a distortion 
in which case internal forces will be developed to counterbalance the external 
forces. If the internal and external forces are balanced, we say that the body is 
in static equilibrium. Otherwise the body can be accelerated which would give 
rise to inertia forces [581. Using d'Alembert's principle, these forces may be 
included as part of the body forces [48J such that the equilibrium equations can 
be satisfied. If the body gets deformed then the deformation can either be 
elastic or not and is subjected to the material properties of the body such as its 
elasticity and its rigidity. If the internal forces induced by the material properties 
of a body are too weak to counterbalance the external forces, then the body 
can be permanently deformed [52]. 

[00267] We assume herein the material to be isotropic with respect to 

some mechanical properties. We then suppose that the material properties are 
the same in all directions for a given point [60]. We also consider an 
homogeneous material which means that its properties are identical at all 
locations. 

The Concept of Stress at a Point 

[00268] Let us consider an isotropic and homogeneous body B. Let 

us assume that B is subjected to arbitrary surface and body forces such that B 
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is in static equilibrium. Let P be a interior point of B and S be a plane surface 
passing through P. S will be referred to as the cutting plane and is defined by 
the unit normal vector n-(ni, n 2 , n 3 ) T . Then S partitions B into two sections I 
and II as shown in Figure 35- 

[00269] Let us assume that AS is a small element of area of the 

cutting plane surrounding P (see Figure 36). 

[00270] Since the body is in static equilibrium then the force system 

acting on each part I and II taken alone must also be in equilibrium. This 
generally requires that some internal forces are transmitted by part I to part II. 
These forces are not necessarily distributed uniformly on every part of the 
cutting plane. Thus they may vary in magnitude and direction over it. We 
generally want to determine precisely that force distribution at every point of 
AS . The term stress is used to define the intensity and the direction of the 
internal force 4f acting at point P. Using the Cauchy stress principle [60] we 
define the stress vector (or traction vector or traction forces) t n at P as 

t tt = lim 

assuming that P remains an interior point of AS as its area reduces to zero. 

[00271] Let us mention that t n is not necessarily in the direction of the 

normal vector n at P. However, it may be decomposed into a component 
perpendicular to the cutting plane, called the normal stress, and a component 
parallel to it, called the shear stress. The normal stress attempts to separate 
(bring closer) the material particles after a compression (an extension) of an 
elastic body when it tries to recover its original state. On the other hand, the 
shear stress acts parallel to the cutting plane and tends to slide adjacent planes 
with respect to each other (see Figure 37). 

[00272] We must notice that the stress vector t° is defined with 

respect to the cutting plane's unit normal vector n. Since there are infinitely 
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many cutting planes going through P there are also as many stress vectors 
defined at P. Juvinall [58] defines the state of stress at a point P as a complete 
description of the stress magnitude and direction for ail possible cutting plane 
passing through P. Fortunately, this description can be fully obtained by 
considering any three mutually perpendicular planes passing at P [58]. For the 
sake of simplicity, we usually use the three axes defined by the three canonical 
vectors Xi, X2 and X3. 

[00273] Let us define <r 9 as the stress component in the direction of xj 

when the normal vector is parallel to the axe defined by xj. If i = j then a„ 
represents a normal stress. Otherwise, a v is a shear stress. With these 
conventions, the component V in the direction of * of the traction force t n 
depend on the normal stress cr u , the shear stresses cr /; and a a and the normal 
vector n such that (see Figure 38) 

tf = <Jun\ + £72»Tl2 + (731^3 
3 

= 

J=i (64) 
Equation 64 is known as the Cauchy stress formula. 

[00274] Since each of the three coordinates axes involves six , there 

is a total of nine stress components. However the equilibrium of moments at P 
[49] gives that only six of these are independent that is <r v = a p for all i, j = 1 , 2, 
3. This means that the state of stress at a point is fully determined by <r„ , <r n , 
<r 3J . o n% er„and tr a . 

The Concept of Strain at a Point 

[00275] Any non rigid body goes through deformations and distortions 

when subjected to forces. The body can either extend or contract (deformation) 
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or have a geometric modification of its shape (distortion). Figure 39 presents 
these concepts. 

[00276] The term strain refers to the direction and intensity of the 

deformation at any given point with respect to a specific plane passing through 
that point [58]. As for stress, the strain is defined according to a specific cutting 
plane. The state of strain is defined by Juvinall [58] as the complete definition of 
the magnitude and direction of the deformation at a given point with respect to 
a all cutting planes passing through that point 

[00277] As for the state of stress, the description can be obtained by 

considering any three mutually perpendicular planes passing at P. We can 
therefore see a great similarity between stress and strain. However, there is a 
major difference between them: strains are generally some directly measurable 
quantities while stresses are not. Fortunately, stresses can be computed from 
strains (and vice versa) using a constitutive equation. 

[00278] As for stress, two types of strains can be defined. First, there 

are the strains which result of a change in the dimensions of the body 
(deformation). Let B be the same body defined hereinabove, AS be a small 
element of B of length Ax, in the direction of x t and AB' be the deformation of 
AB such that Am, is the change in lenght of AB after the application of a force 
in the direction of Xj (see Figure 40).The normal strain e a at P in the Xi direction 
with respect to a cutting plane having xj as normal vector is the unit 
deformation of a line element [65] in the direction of X|. It is formally defined as 

Aui dm 

€u — lim -r = 

[00279] The normal strain is clearly the unit change in length per 

unit original length for the element in the direction of x». Since it is the ratio of 
two units of length, it is dimensionless even if it is sometimes expressed as 
units of length per unit of length such as inches per inch. 
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[00280] Let us now suppose that we have two perpendicular lines PB 

and PA of length Ax 7 and Ax* respectively in the direction of xj and x k (see 

Figure 41). 

[00281] Let us assume that after a distortion points A and B move 

respectively to A 9 and B' while P remains fixed. The lines PA and PB have 
been rotated of angles 0 A and 0 W such that 

^ = tan(^ fc ) and ^ = tBn(0 kj ) 
Ax k 

where Aw, and Aw* are respectively the displacements of B and A in the xj and 
xk directions. 

[00282] If we assume that only small distortions occur, then we can 

approximate both tangents by their angles. Thus 

0 jk tan(^) = and 0 kj ^ tan(O kj ) - ^ 
or, taking their infinitesimal analogous 

e jk = hm "t = and 0 kj = lim -r— = -x— 

[002B3] The shear strain at P with respect to the cutting plane 

having xi as normal vector is the angle in radians through which two orthogonal 
lines in the undistorted body are rotated by a distortion [56]. That is 
7ik =0jh The toto subscripts in y ft have a similar meaning as for stress. 
For instance, y & is the strain acting on two adjacent planes perpendicular to the 
xj axis and sliding them relative to each other in the x k direction. 

[00284] Let us recall that these definitions have been made under the 

assumption that only very small displacements occur in the body. The normal 
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and shear strains are supposed smalt compared to unity [50]. If we relax this 
constraint in order to include large deformations then the system to solve for 
the computation of the forces, the stresses, the strains or the displacements 
becomes non-linear and then harder to solve. This is sometimes necessaty in 
some problems where large deformations can occur such as for thin flexible 
bodies [50] of for the moderation of human tissue [53]. However, if we restrict 
ourselves to small deformations, then the approximations made to define the 
shear strains do not induce too much errors and are widely accepted in the 
classical theory of elasticity [69, 65, 49, 50]. 

[00285] Finally, Equation 65 clearly shows that the shear strain is also 

a dimensionless quantity since it is the ratio of two units of length. If we define 
for i * j 

e«« \lH 1,2,3) 

then we get the strain-displacement relations (or the kinematical relationship 
[69]) 




[00286] As for stresses, there is a total of nine strain components at 

each point of the body (three per mutually perpendicular cutting planes) but by 
symmetry they can be reduced to six, that is m Yv for all j, k = 1 , 2, 3 with j * 
k. The state of strain at any point can then be described by e u , e n , e S3 , e n , 
£, 3 and e n . 

Relations Between Forces, Stresses and Strains 

[00287] As mentioned hereinabove, strains are measurable quantities 

while stress are not even if both can be computed from the other. This is due to 
the fact that some knowledge about the material of a body is necessary to 
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measure the stresses from strains and vice versa. For instance, a steel beam 
and a rubber beam induce different internal forces when bent equally. 

[00288] We will now fill this gap between strains and stresses by 

stating the physical laws relating them. This hole between strains and stresses 
need to be filled by a constitutive equation (or material law) which reflect the 
internal constitution of the materials. The material law for the linear elasticity 
problem is known as the Hooke's law. 

[00289] Before stating the law, let us recall that a material has an 

elastic behavior when it satisfies the two following conditions: 

1 . The stresses depend only on the strains. 

2. Its properties allow a body to recover its original shape when the external 
forces applied on the body are removed [60]. 

[00290] If theses conditions are not satisfied, a body is said to have 

an inelastic behavior. Any body may be seen as having an elastic behavior as 
long as it is not deformed beyond a limiting value. This value is called the 
elastic limit [52, 49] and is usually defined as the maximum value of stress that 
a body can undergo without permanently being deformed. 

[00291] In addition, if the stress is a linear function of the strain we 

say that an elastic material has a linear elastic behavior. In what follows, we 
assume that the material has a linear elastic behavior. 

[00292] As mentioned in [60] and [49], in many situations the problem 

of elasticity can be considered as a 2D problem. This particular rase is known 
as plane elasticity. Two basic types of problems compose the plane elasticity. 
The problems in which the stress components in one direction for a body are all 
zero are refered to as plane stress problems. On the other hand, if all the strain 
components in one direction for a body are zero then the state of strain for that 
body is refered to as plane strain (see [65, 60, 58]). 
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[00293] We consider our problem as a plane strain problem. This 

distinction is important since the constitutive equation slightly differs according 
to the fact that we consider a plane stress or a plane strain problem. 

The Hooke's Law 

[00294] When a rubber ball is compressed its diameter in the 

directions perpendicular to the applied force gets larger. A similar phenomenon 
happen occur when a rubber band is extended and its cross section gets 
smaller 1 . In fact, these changes in dimensions happen in all materials even if 
they cant always be noticed by a naked eye [52]. 

[0029S] When a stress is acting on an isotropic and homogeneous 

body in only one direction (uniaxial stress), one can show that the transverse 
strain s x (the strain in a perpendicular direction) is directly proportional to the 
strain s induced by the stress 

e 1 - = —ve 

The ratio 




is called the Poisson ratio. It is supposed constant when the stress is below the 
elastic limit [66]. 

[00296] The linear relationship between a uniaxial stress cr H in a 

direction X| and the corresponding strain e a is known as the Hooke's law [52, 
58] and is written as 




1 Example taken in [52] 
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where E is the Young's modulus of elasticity. Of course, the Hooke's law is 
valid if the stress is not beyond the elastic limit of the material. 

[00297] As pointed out by Boresi [49], the stresses at any point 

depend on all the strains in the neighborhood of that point. Thus the total 
deformation in the xi direction depends not only on the stress in that direction 
but also of the deformations in the other two perpendicular directions. For 

instance the normal strain s fi does not only depend on ^(Hooke's law) but 

also of the transverse strains s n and % such that the total deformation in the 
direction of X\ is 

- a * „?M v akk 
E E E 

= ^;[(Tii-v(<rjj + akk)] (i^J^k) 
& (67) 

[00298] Equation 67 is the normal strain-stress relation. For isotropic 

materials, it can be shown that the normal strains are not influenced by the 

shear stresses [52]. Consequently, the shear stresses only induce shear strains 

and they are related by the relation 

G x * JJ (68) 

where G= — - — is called the modulus of rigidity. Equations 67 and 68 are 
2(1 +v) 

known as the Generalized Hooke's law for linear elastic isotropic materials [52]. 
These equations can be inverted in order to express the stresses as functions 
of the strains 
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E 



°« " ( i + y) (l-2 y ) [(1 ^^ + y( ^ +g * fc)1 

(69) 

= *»«*-if?* (Mi) (70) 

Relation Between Forces and Stresses. 

[00299] It is well known that the conservation (balance, equilibrium) 

laws constitute an important class of equations in continuum mechanics. They 
relate the change in total amount of a physical quantity inside a body with the 
amount of this quantity which flows through its boundary. These laws must be 
satisfied for every continuous materials. Local differential equations are usually 
used to express these laws. In what follows, we present the linear momentum 
which is relevant for the linear elasticity problem. 

[00300] To every material body B is associated a measure of its 

inertia called the mass. This measure may vary in space and time inside a 
body. Let V be the volume of B, S its bounding surface and Am be the mass of 
a small amount of volume 6V . We call mass density 

p = p(x,t) = lim ~j 
r r\ > * av-x) AV 

[00301] Let us assume that distributed body forces pb t and tractions 

forces ti n are applied to S (see Figure 42). Let also assume that B is moving 
under the velocity field vj = vtfx.t). The quantity 



P*(t) = /// P»idV 
v 



is called the linear momentum of B. The principle of linear momentum [49] 
states that the resultant force acting on a body is equal to the time rate of 
change of the linear momentum. Thus 
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S V 



Forces acting on the body (7.^ 

[00302] Recalling Equation 64 



3 



where n = (ni» n 2 , ri3) T is the unit normal vector to the surface and 
a = (<r u ,<r„,<rj,) r and using Gauss's divergence theorem, we have 



eft 



/// < M < dV = Iff %l"-ff J 

V v V (72) 

[00303] Since V is arbitrary then the integral sign can be retrieved 

leading to the local equations of motion 

p— = V -a + pbi 
dt (73) 

[00304] The global equilibrium equations can be obtained assuming a 

zero velocity field in Equation 72 

//TV* dV + ////*< dV = 0 (i = 1,2,3) 

v . v 



> „ ^ 



Internal forces External forces (74) 

and their local counterparts are 

V-<r + pD<=0 (2 = 1,2,3) (?5) 



— m 
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[00305] Let us now summarize the decomposition of the problem. We 

introduced the relation between the global displacement U of a body and the 
corresponding strains (see Figure 43(a)) hereinabove. We also presented the 
constitutive equation relating the strains and the stresses (see Figure 43(b)). 
Finally, we described how the stresses are related to forces using the linear 
momentum principle (see Figure 43(c)). A general scheme similar to the one 
presented by Tonti [72] may then be introduced to summarize how the internal 
reaction forces of a body are related to the global displacements of that body 
(Figure 44). 

Discrete Representation of Images 

[00306] Some algebraic tools used to model the image will now be 

recalled from the above description. An image is composed of two distinctive 
parts: the image support (pixels) and some field quantity associated to each 
pixel. This quantity can be scalar (e.g. gray level), vectorial (e.g. color, 
multispectral, optical flow) or tensorial (e.g. Hessian). We model the image 
support in terms of cubical complexes, chains and boundaries. With these 
concepts, we are able to give a formal description of an image support of any 
dimension. For quantities, the concept of cochains which are representations of 
fields over a cubical complex is introduced. For the use of these concepts in 
image processing, see [45]. 

[00307] An image is a complex of unit cubes usually called pixels. A 

pixel y c W is a product 

7 = / l xJ 2 x,.,x/ n 

where lj is either a singleton or an interval of unit length with integer endpoints. 
Then lj is either the singleton {k} and is said to be a degenerate interval, or the 
closed interval [k, k+1] for some *eZ. The number ?e{0,l,...,/i}of non- 
degenerate intervals is by definition the dimension of y which is called a q- 
pixel. Figure 45 illustrates three elementary pixels in ft 2 . 



• 
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[00308] For q > 1 f let J « {k^k^k^} be the ordered subset {1, 2, 

n} of indices for which I kj =[a J9 bj] is non-degenerate. Define 

A ki a = h x . . . Jj^i x {a,} x I kj + X x . . . x I n 

and 

B^a = h x . . . J^-i x {fy} x J fcj+1 x ... x J n 

[00309] The y^and the B k/ are called the (q-l)-faces of <r. One can 

define the (q-2)-faces, .... down to the 0-faces of a the same way. The faces 
of y different from y itself are called its proper faces, 

[00310] By definition, a natural orientation of the cube is assumed for 

each pixel. Suppose that y denotes a particular positive oriented q-pixel. It is 
natural to denote the same pixel with opposite orientation by -y . Examples of 
orientations are given in Figure 45. A cubical complex in <R" is a finite collection 
K of q-pixels such that every face of any pixel of the image support called K is 
also a pixel in K and the intersection of any two pixels of K is either empty or a 
face of each of them. For example, traditional 2D image models was only 
considering pixel as a 2D squared element. Definitions presented before allows 
us to consider 2-pixels (squared elements), 1 -pixels (line elements) and 0- 
pixels (punctual elements) simultaneously. 

[00311] In order to write the image support in algebraic form, we 

introduce the concept of chains. Any set of oriented q-pixels of a cubical 
complex can be written in algebraic form by attributing them the coefficient 0,1 
or -1, if they are not in the set or if they should or not be taken with positive 
orientation, respectively. In order to represent weighted domains, we must 
allow arbitrary integer multiplicity for each q-pixel. 

[00312] Given a topological space X in terms of a cubical 

complex, we get a free abelian group C q (X) generated by all the q-pixels. The 
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elements of this group are called q-chains and they are formal linear 
combinations of q-pixels [45]. A formal expression for a q-chain Cq is 

c 9 =Hr tmK *iri where ^eZ. 

[00313] The last step needed for the description of the image plan is 

the introduction of the concept of boundary of a chain. Given a q-pixel y , we 
define its boundary dy as the (q-1)-chain corresponding to the alternating sum 
of its (q-1)-faces. The sum is taken according to the orientation of the (q-1> 
faces with respect to the orientation of the q-pixel. We say that a (q-1)-face of 
y is relatively positively oriented with respect to the orientation of y if its 
orientation is compatible with the orientation of y . By linearity, the extension of 
the definition of boundary to arbitrary q-chains is expedient. For instance, in 
Figures 45(b) and 45(c), the boundary of the 1 -pixel a is x 2 - xi and the 
boundary of the 2-pixel A is a+b-c-d and we say that a and b are positively 
oriented with respect to orientation of A but c and d are negatively oriented with 
respect to orientation of A. Let us notice that the boundary of a 1 -pixel is always 
the difference of its boundary points. The boundary can be defined recursively. 
Suppose a (q-l)-chain and a q-chain y g defined as r« -j^,x[a,6], the 

boundary of y q can be recursively written as 

dy„ = ^r H x fc b] + (-l^Ofc-i x {6} - 7 ,_ t x {a}) (?6) 

[00314] In order to model the pixels quantity over the image plane, we 

are looking for an application F which associates a global quantity with ail q- 
pixels y of a cubical complex. We note this <F,y>. This quantity may be any 
mathematical entities such as scalars, vectors, etc. For two adjacent q-pixeis /, 
and y 2 , F must satisfy < F t + Iff* >= Aj <F % y x >+X l < F, y 2 > which means 
that the sum of the quantity over each pixel is equal to the quantity over the two 
pixels. The resulting transformation F\C q (X)-± « is called a q-cochain and is 

used as a representation of a quantity over the cubical complex. 
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[00315] We finally need an operator which associates a global 

quantity to the (q+1)-pixels according to the global quantities given on their q- 
faces. Given a q-cochain F, we define an operator 5, called the coboundary 
operator, which transforms F into a (q+1)-cochain dF such that 

< ST, 7 >=< Bj > (77) 
for all (q+1)-chains y. The coboundary is defined as the signed sum of the 
physical quantities associated with the q-faces of y. The sum is taken 
according to the relative orientation of the q-faces of the (q+1)-pixels of y with 
respect to their orientation. Figure 46 presents an example of the coboundary 
operation for a 2-pixel. 

Representation of the Equilibrium Equation 

[00316] The basic laws of Figure 44 with concepts of algebraic 

topology in order to get a generic algorithm for solving the equilibrium equation 
74 will now be modeled. 

[00317] The algorithm is resumed as follows: 1) The image support is 

firstly subdivided into cubical complexes; 2) Global quantities are computed 
over pixels of various dimensions via cochains according to basic laws; 3) The 
constitutive equations 69 and 70 are expressed as a linear transformation 
between two cochains. 

The Relative' Displacement 

[00318] Let B be a body in a 3D space and K p be a 3-complex 

representing the subdivided spatial support of B. Let us consider a 0-cochain V 
and a 1-cochain D such that D is the coboundary of U 



7 h-> <X>,7 >=< <5W,7 >=<U,dy > 
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[00319] Figure 47 presents some examples of V and D for a 3-pixel of 

K", 

[00320] We must now specify the computational rules for both 

cochalns V and D. Recalling the strain-displacement relation (Equation 66) 



* y 2 [dxi ^ dxjl 



(79) 



we have an application s' 



[00321] 

define 



[00322] 



U £*(U) = (£u,S22 } C33>£l2,ei3»^23) T 

Omitting the shear strain components as in [67], we may 
e : R 3 -* 

U »-+ e(U) = (fn,€:22,£33) T = VU (8Q) 
Using the global form of Equation 80 over a 1 -pixel y D such 



that dy D = x, - x„ , we have 

/ e(U)d7D= r*e(U)rf7i>= f VUd7D 
/*> ^# J *# (81) 

where is an infinitesimal part of the domain y D . Since Vu is a 

conservative field, we can apply the line integral theorem [55, 71] which states 

that for a conservative field F(x) = Vf(x) and for two points A and B in an 

opened connected region containing F(x) the integral of the tangential part of 

F(x) along a curve R joining A and B is independent of the path 



[00323] 



/"*F(x)di2 = /(B)-/(A) 

J A 

From Equation 81 , we then have 
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f 



VUd7i> = U(a;,)-U(x # ) 

(82) 

[00324] On the other hand, applying the cochain D to the 1 -pixel y D 

we also have 

< 2>, Id >=< U, dj D >= U(x+ - x # ) = U(x.) - U{x#) (fl3) 

which is the same form as Equation 82. We then define U(x) - U(x) 
Consequently, the location of the displacement vector U must correspond to 
the 0-pixels of K" . The previous definitions are extended by linearity to the 1- 
chainsof K p '. 

The Force-Stress Relation 

[00325] Let us consider another 3-complex J?* also representing the 

subdivided spatial support of the body B. Let us also consider a 3-cochain 3 
and a 2-cochain S such that.? is the coboundary of S 

F:C3(K?) -» K 3 

7 -> <^7>=<^,7>=<^,^7> (84) 

[00326] Figure 48 presents some examples of 3 and S for a 3-pixel of 

K", 

[00327] Let y F be a 3-pixel of K'and y s be a 2-chain over K' such 

that Ys = QYf • Let us assume that the 2-faces y Sj of y F are relatively positively 
oriented with respect to the orientation of y F . By the definition of the 
coboundary we have 



ft 



(85) 
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[00328] Again, we want to determine the computational rules 

associated with F and S. If we set a zero velocity field in Equation 71 , we have 



II ^^ r== If ~ a * ' ndS 



v s 
where cr t = (cr u ,<x 2 .,<7 3f ) . To fulfill Equation 85 we define 



<^,7f >= fff pbidV 



(86) 



and 



<S iy j Sj >- ^ - ndS (i= 1,2,3) 

S (87) 

where 
and 

The Stress-Strain Relation 

[00329] We have presented exact global versions of Equations 66 on 

K p and 74 on K'by the means of Equations 83, 86 and 87. In order to 
complete the scheme of Figure 44, we need to represent the Hooke's law 
(Equations 69 and 70) which links the local values of *(x) and a(x). Since 
Equations 69 and 70 are constitutive equations, we can't provide a topological 
expression of them. Instead we express them as linear transformations 
between the cochains D and S 
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[00330] To find this transformation, let us recall Equation 87 

< Si,-Ysj>= J J ^-<Ti*ndS (* = 1,2 ? 3) 

which links the cochain S with the strains s using the generalized Hooke's law. 
Unfortunately, the strains are only known at a finite number of points and must 
then be approximated over the whole domain S with an approximation function 
*(x). We choose e(x) such that for each 1-face y D of a 2-pixel y of K p w © 
have 



e(x) • dR=< V^o > 
*"> (88) 

where dR is an infinitesimal part of the domain represented by y D . Let us 

notice that we only need to approximate the normal components of e . In fact, if 

we have 

VU(x) = (f£(x), f^(x),f|(x)) = (eu(x),£ 22 (x) J ? 33 (x)) T = ?(x) 

where w(x) is the approximated displacement vector over y and if s(x) is 

chosen to satisfy Equation 88, then the vector t7(x) is fully determined. Then 
the shear components of e(x) can be computed by appropriately differentiating 
the components of C/(x). Using this remark and applying the generalized 
Hooke's law to e (x) satisfying Equation 88, we have 

E 

= (1 + y)(l - 2u) ^ ~ + y <5i( x ) +g fcfc(x))] 

E 

= (l + vf vW = 1 > 2 ,3) 

at all point of y . Equation 87 is then replaced by 
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< Si>7Sj >= // -a,(x) ■ ndS = A*(e) (i = 1,2,3) 

JJ 8 (89) 
which depends on the choice of the approximation function e(x) and of the 
relative position of K s with respect to K p . 

Summary of the Algorithm 

[00331] We now summarize the algorithm used to find an expression 

of the internal forces according to the displacements of a body. The input data 
for this algorithm are the cochain U and the material properties of the body (the 
values of E and v). 

1 . Choice of the location of K p with respect to K* . 

2. Computation of the cochain D. 

3. Computation of the cochain S 

(a) Choice of the approximation function e (x). 

(b) Application of Equation 89 to express S as a function of the 
displacement components. 

4. Computation of the force by applying Equation 85. 

Applications. 
Active Contours 

[00332] We apply the above described approach to a 2D active 

contour model based upon a Lagrangian evolution of the curve S 

|? + KS = <?) 

^ (90) 
where K is the matrix which contains the regularization forces of the curve. 
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[00333] This dynamic system is discretized in time using a finite 

difference scheme. For a given time step Af , we can approximate the time 
derivative by 

OS ^ St+&t — St 
dt ~ At 

[00334] The curve deformation is governed by each vertex 

displacements compared to their neighbors until the equilibrium between the 
inertia forces, the image forces and the internal forces. We solved Equation 90 
using an explicit scheme 

St+Ai = S t + At(F„ t - KS t ) (g1) 

[00335] Assuming that the initial curve So was in an equilibrium state 

and that the initial body forces F 0 = KS 0 are constant during the deformation 
process, we can add these forces to the external forces F ex t leading to a 
modified version of Equation 91 

= S t + At(F„ t + F 0 - KSt) 

= S t + At(F m - (KS t - KS 0 )) 

= S + AttF^-KU) (92) 

where U is the displacement vector of the curve S. 

[00336] The image subdivision process is similar to the one 

presented in [57]. We want to solve Equation 92 for local U(x) located at the 
center of each pixel and known initial curve So closed to the solution. Following 
the steps presented hereinabove, we first position the two dimensional cubical 
complexes K p and K* . As mentioned, K p must be placed in order to have its 
0-pixels corresponding to the center of the image pixels. We positioned K* in 
such a way that its 2-pixels coincide with the image pixels. Thus, the 2-pixels of 
AT'are rectangular and symmetrically staggered with the 1 -pixels of K p and 
each 1 -pixel of intersects orthogonally in the middle of a 1 -pixel of K p . 



— • • 
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Mattiussi [61] showed that this way of positioning fallows the use of lower 
order approximation polynomials without losing accuracy. Figure 49 shows the 
two complexes positions for a 5x5 image. 

[00337] In order to solve Equation 92, we need global values F over 

each pixel of K s . Since these values are generally known, we did not try to 
express them in a topological way. In our examples, we use the gradient field of 
the bright line plausibility image obtained using a line detector proposed in [54]. 
We assume that the gradient provides global values valid over the whole pixel. 
Thus we set F 3 VL where L is the line plausibility image, VL = g'^L and g' a 
is the Gaussian derivative at scale a. We also choose an approximation 
function s (x)=(s u (x), £b(x)) t . For simplicity, we assume that (x) = V£/(x) 
arises from a bilinear approximation 

U(x) = (Ui (x), U 2 (x)) T = a + hx x + cx 2 + dx x x 2 
[00338] Thus we have 

fu(x) = b + dap 2 
£ 22 (x) = c + dxj 

[00339] Since s u (x) and ^(x) must satisfy Equation 88 for all 1- 

faces y D of any 2-pixel y of K p as in Figure 50, the following relations must 
hold 




T>2 = / *(A,a: 2 ) • j dx 2 

T>i = / e(x u A) ■ i dxi 
Jo 
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V 4 = /*e(0,A) • T<**2 
Jo 



from which we obtain 



[00340] From Equation 93 and the definition of the normal strains, it is 

straightforward that 

U(x) = k + i (T> x x t + X> 4 tf 2 ) + i (I>3 - A + Z> 2 - ^4) *ia; 2 

A A z (94) 

where k=U(0) . Using Equations 83 and 94 we have 

U(x) s U(*i,* 2 ) = U(0,0) + ^(U(A,0)-U(0,0))xi + i(U(0,A)-U(0,0))x 2 
< U (°> °) + U < A > A > ~ U (°> A > - U < A > °)) 

(95) 

from which we can deduce the values of 07 (x) = ( <x w (x), o 2i (x)) T . 

[00341] The last step is the computation of the internal forces F for 

each 2-pixel of K* . With K p andJT positioned as mentioned before, we have 
that each 2-pixel y F of K 9 intersects four 2-pixels y A% y B , p c and y D of K p . 

That is, we must consider four approximation functions of , of , of and of 

corresponding to the four intersecting 2-pixels of K p (see Figure 51) 

[00342] We find the value of the cochain S over the four 1-face of 

yby the appropriate integration 

Sf = j*of --T*i + /° A 3f --t^ 
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s? = 



[00343] 



Using Equation 85 we have 



< Fi , 7F > = S\ + 5? + Sf + 5/ 



(96) 

[00344] Substituting Equations 96 and 95 in Equation 96, we can 

express the internal forces F as a function of the displacement U. As an 
example, we present the values of F for the 2-pixel y F of Figure 50 with A = 1 

F x = C{(3 - 4i/)u-i.i + (2 - 8^)^0,1 + (3 - 4f/)u la + (10 - 8f)u_i l0 + (-36 + 4%v)u 0 $ 
+(10 - 8i/)ui, 0 + (3 - 4i/)u-i.-i + (2 - 8i/)u0, -1 + (3 - 4u)ui^i - 2v- U 
+2t?, a -f 2t/_ x ,_ l -2« li - 1 ] 

F 2 = C [(3 - 4i/)u_i,i + (10 - 8i/)«o,i + (3 - 4i/)u ltl + (2 - 8i/)u_ li0 + (-36 + 48i/)u D # 
+(2 - 8i/)m >0 + (3 - 4^)w-i,-i -f (10 - 8i/)u 0 ,-i + (3 - Iv - 2v_ M 
+2v M + 2t/_,,_ t - 2V!,-,] 

where 



_ b 

^ ~ 16(l+i/)(l-2i/) 

and 



"-[:] 



[00345] Equation 97 induces a linear relationship between a pixel and 

its neighbors. This relation is used to build the stiffness matrix of Equation 91. If 
we let (i=1, 2) be the displacement vector for the component X| then we 

have 
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JV *2 



E 



[00346] 

stiffness kernels. 



16(l + i/)(l- 


-2u) 


E 




16(l + i/)(l~ 


2u) 


E 




16(l + */)(l- 


2v) 


E 





3- 


-4z/ 


2-8i/ 


3 


-Av 


10 


-81/ - 


36 + 48^ 


10 




3 


- Av 


2-8i/ 


3 


-Au 


-2 


0 2 " 








0 


0 0 








2 


0 -2 








-2 


0 2 








0 


0 0 








2 


0 -2 









3-4*/ 

2- 8i/ 

3- 4r/ 



10 - 8u Z-Av 
-36 + 481/ 2 - 8i/ 
10-8i/ 3 - Au 



16(l + i/)(l-2i/) 
The pairs (iv;,^), (i = 1, 2) will be referred to as the 



Computation of the Displacement Vector 

[00347] The assumption made when calculating the displacement 

vector will now be explained. Let v be a vertex of s subdivided curve S and v' 
be its corresponding vertex in the deformed curve S' . Let us denote by U[v] the 
entry in the displacement vector U corresponding to v. Let us suppose that the 
displacement is constant in each direction everywhere that is U[v] = (ki, k 2 ) T 
with ki, k2e SR for all v in S. Since the sum of all entries of either N 1 . JV' N 2 
or A£is zero, it follows that Fi[v] = F 2 [v] = 0 for all v in S which means that 
there is no internal force induced. The computation of the internal forces with 
the stiffness kernels is then invariant with respect to translation. 



[00348] Let v 1f v 2 , v 3 , v 4 and v 5 be five adjacent vertices of S and v[ , 

v 2 . V 3 , v' A , v' s be their corresponding vertices in S' (see Figure 52). 



• 
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[00349] Let v ifXj be the xj coordinate of the spatial representation of 

the vertex V|. Then we have 

[00350] By the translation invariance property, we have 

where [v 2tXi -v 2 ^] stands for a matrix whose all entries equal v' 2tXx -v 2tXr The 

displacement component used to compute the internal force Fj at vertex v 3 is 
then 

which is the relative displacement of the vertex v 3 with respect to v 2 . However, 
nothing prevents us from computing the relative displacement of v 3 with respect 
to v 3 . In this case, the x k displacement component used to compute the internal 
force F| at vertex v 3 would be 

u* k M = («^-^*)-( v S«i-»4») 

[00351] In order to take into account these facts, we use the average 

value of the relative displacements for all adjacent vertices to V2. Thus 

U* k [v 2 ] = |[K, % -tkJ-fo^ 

(97) 

[00352] Reorganizing the terms of Equation 97, we have for an 

arbitrary vertex v s of S 
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U Xk [vi] = ^+^- 2 ^+^-i^ j _ ^i+l* k ~ K* t ± "Li** 



or 

\j[ Vi ] = - 2 ^' + t, »-i j _ -2^H-^_ t j 

2[dv* dv*\ (98) 

if we assume a finite difference approximation of the second derivative of S and 
5" with respect to their vertices. 

[00353] Let us finally notice that the second derivative of the curve S 

in Equation 98 can be computed using Gaussian derivatives 

where a controls the degree of smoothing. Such a computation of the second 
derivative of S allows to obtain smooth results by simulating a smoother target 
curve. 

Experimental Results 
Active Contours 

[00354] The approach proposed herein has been experimented on 

real and synthetic images in the context of high-resolution images of road 
databases. For each image, we have compared our results with another 
method. 

[00355] Figure 55(a) presents the results for our physics-based 

method (PB) for an aerial image while Figure 55(b) shows the results for the 
finite element (FEM) method (a=??? f p=???). We simulated a material similar 
to rubber (see the Table in Figure 53) with E=150 and v=0.45. In both images, 
the image force is the gradient (<r=1.5) of the bright line plausibility image 
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obtained using a line detector proposed in [54] with the line detection scale set 
to 0.8. Figures 54(a) and 54(b) show respectively the initial curve So and the 
bright line plausibility image. 

100356! Figure 56 shows a SAR image in which the initial curves are 

drawn in white. The line plausibility image (<y=1.5 and line detection scale 
=0.8) of both bright and dark lines is shown in Figure 57. Figures 58 and 59 
presents respectively the results obtained with PB (E^iso and v=0.45) and 
FEM (a=??? and p=???) methods. One can notice that the PB curves are 
closer to the shore than the FEM especially in region of high curvature. 

([00357] Figure 60 shows some initialization for the first band of a 

Landsat 7 image (courtesy of XXX). The bright line plausibility image (a =1.5 
and line detection scale =0.8) is shown in Figure 61. Figures 62 and 63 present 
the results obtained with the PB and FEM methods respectively. 

[00358] We have discussed the fact that the deformations obtained 

using our model have a physical interpretation. We also discussed the fact that 
the objects modeled using the PB method have their own physical properties 
and the ability to recover their original shape when the external forces applied 
on them are removed. To illustrate this fact, Figures 64(a) and 64(b) show 
some initialization and the corrected curve for a synthetic image. Figures 65(a) 
through Figures 65(d) show the evolution of this corrected curve when the 
external forces are removed. Figure 65(d) presents both the final curve (in 
black) and the initial curve (in white). One can clearly notice that the curve has 
recovered its original shape but has also experienced a spatial shift. 

Conclusion 

£00359] A new model for active contours was presented. The 

proposed approach decompose the image using an image model based on 
algebraic topology. This model uses generic mathematical tools which can be 
applied to solve other problems such as linear and non-linear diffusion and 
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optical flow [57], Moreover, our model works with either 2D or 3D images and 
can easily be extended to active surfaces and active volumes. 

[00360] Our approach presents the following differences with the 

other methods: 1) We use both global and local quantities; 2) Our model is 
based upon basic laws of physics. This allows us to give a physical explanation 
to the deformation steps; 3) The curves and surfaces have physical behaviors 
such as the ability to recover their original shape once the applied forces are 
removed; 4) We make approximation only when the constitutive equation is 
involved. 

[00361] Although the present invention has been described 

hereinabove by way of preferred embodiments thereof, it can be modified, 
without departing from the spirit and nature of the subject invention. 
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